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L INTRODUCTION. 



This report presents some navigation and almanac algorithms with implementation pro- 
grams that are written in the Radio Shack TRS-so MODEL 4 BASIC language. This language 
can be customised for the Navy adopted Sharp PC-lSOOA or other BASIC language com- 
puters. 

Because the PC-I600A has limited memory, many of the classical navigation formulas 
from spherical trigonometry have been rewritten in a form which minimises or eliminates 
the need to determine special cases — especially those of quadrant. Details are contained 
in Section IL 

The navigation algorithms are detailed in the BASIC program NAVALGOR, presented in 
Section IlL Programs for the determination of the position of the Sun, Moon, Venus, Mars, 
Jupiter and Saturn are contained in the BASIC program NAVEPHM which is presented in 
Section IV. These programs reproduce the tables in The Nautical Almanac and The Air 
Almanac to within 0.2'. 

n. THE GENERAL SPHERICAL TRIANGLE 
AND COORDINATE TRANSFORMATIONS. 

A. Background. Many of the algorithms and procedures used for navigational compu- 
tations were designed for use with tables of sines, cosines and tangents. To save space and 
eliminate repetition, these tables usually contained only values for the first quadrant. To 
generate values for other quadrants, many cumbersome rules were adopted. These rules — 
contained in classical references such as Bowditch [Ref. 1|— are promulgated in modem 
computer programs, causing needlessly additional programming effort. This additional ef- 
fort can largely be eliminated by using the concept of the general spherical triangle which 
is described in the next section. 

B. Theory. The formulas for the general spherical triangle expounded on by Wm. Chau- 
venet [Ref. 3| in his book, A Tnatice on Plane and Spherical Trigonometry, which was first 
published in 1850. One usually thinks of a spherical triangle as looking something like the 
object depicted in Figure 1, where the arc length of each leg is less than 180^ Chauvenet^s 
general triangle does not have the 180® limitation and so objects similar to the one shown 
in Figure 2 are also considered to be spherical tiian^es. In Figure 2, the crosspoint P 
is 180® ^m the vertex of angle A In labeling spherical triangles it is important that an 
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Figure 1. A SjAerkal TMangie 




Pigue 3. A General Spherical TMangle 

aa^e be labeled with a capital letter and that the opposite side be labeled with the same, 
but lower case letter. 

Depending upon the quantities to be determined, either the set of equations 

cosa cos 6 cos c + sin 6 sine cos A, 

sin a cos ^ = cos 6 sin c -sin 6 cos c cos A, (1) 

sinasin j? s sin6sin A, 

or the polar form of Equs. (1) 

cosA s — cos^cosC + sin^sinCcosa, 
sinAcos6 = cos^sinC + sin^cosCcosa, (2) 

sin Asin6 = sin^sina. 
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may be used. 

Chauvenet stated that, *It will be foimd that Mike six coses of the general triangle 
admit two sohttiona, bat that thef all become determinate, when, in addition to the other 
data, the sign of the eine or cosine of one of the parts is given,* 

Fonnnlas for spherical coordinates should be developed so that quantities snch as 
latitnde, declination and altitnde, which have values in the range -90^ to +90^, are de- 
tennined by arcsine or arctangent fonnnlas; so that quantities such as longitude, asimuth, 
Greenwich hour angle and local hour angle, which have values in the range —180® to +180® 
or 0® to 360®, are determined by both a sine and a cosine fonnula used in conjunction with 
a quadrant determining arctangent function (the qatn function); and so that a quantity 
such as distance, which has a spherical arc value in the range of 0® to 180® (assuming the 
user wishes the shortest distance), is determined by an arccosine formula. 

Uhutrutiofn 1. Consider the ^inverse* problem of spherical geometry, which is: Given 
the latitude, and longitude. A, of two points Pi(^i,Ai) and i^(^,A 3 ), determine the 
distance d, the forward asimuth ori 3 and the backward asinmth or^i. The solutions will be 
derived using the convention that eastern longitudes are positive and western longitudes 
are negative. Also, northern latitudes are positive and southern latitudes are negative. 
The geometry is illustrated in Figure 3. 



North Pole 




Pignn S. Geometry of the Nari^tion TMan^le. 
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Solve for d and ais. Label the angles and legp so that A = A3 — Ai, B » <213, 
C = 360* - (>31, a = d, 6 = 90*-^ and e = 90* - ^1. Substituting these into Equs. (1) 
and simplifying the differences of angles, we obt^ 

cosd = 8in^8in^i +cos^cos^icos(A3 - Ai), 
sm<icoeai3 = sin^cos^i -C08^sin^iC0s(A3 - Ai), 

8md8inax3 =co8^sin(A3 - Ai). 

The first equation can be used to solve for d. Using the principle angle of the arccosine, d 
lies in the range 0* to 180*. If the non-principle angle is used, d lies in the range 180* to 
360*, we can either travel the shortest great circle distance horn Pi to P2 or we can travel 
the great circle route which goes around the backside of the earth — the former solution is 
assumed to be preferred. Hence, d is restricted to lie between 0* and 180*. Multipfy d by 
60 to get the distance in nautical miles. 

Chauvenet states that this system of equations has two solutions, but has only one 
solution if the sign of the sine or cosine of one of the parts is known. Since d lies between 
0* and 180*, the sign of rind is known, thus the system has only one solution. !\irther, 
and most important, the rign of sin d is positive for all d between 0* and 180* and so there 
are no special cases. 

The first equation is used to compute d: 

d = arccos{8in^rin^i + co 8 ^cos^iC 0 s(A 3 >-Ai)|. (3) 

The second and third equations are used to compute ori3. A natural tendency would be to 
solve either the second equation or the third equation, but doing so loses the information 
concerning the quadrant of 013. Use the second equation to determine cos ori3 and the third 
equation to determine sin 013. With both the sine and the cosine known, the quadrant is 
uniquely determined using the quadrant determining arctangent function. The solution is 

ai3 = qatn(sin ai3, cos ais) 

as long as d is not 0* or 360®. Even though the chance of division by zero with the 
occurrence of a value for d of exactly 0* or exactly 360* is very remote, it can be eliminated 
entirely. Since the sign of sin d is positive or zero in both the second and third equations, 
it effectively cancels out using the qatn function. The asimuth 013 is found solving 

tti3 =qatn{co8«^3sin(A3 - Ai),rin^3C08<^i -coe^sin^i cos(A3 - Ai)], (4) 
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To find tlie back aiinmtli, label the angles and legs so that A = Ai , J9 — 360® -a^i , 
C = ai 3 , a = dy 6 = 90*^ and e = 90® -^ 3 . Substitute these into Equs. ( 1 ) and simplify 
to obtain 

asi = qatn[-co 8 ^i sin(Ai - A 3 ), sin cos ^ - cos sin ^ cos(Ai - A 3 )]. (5) 

ISastratioii 2 * Consider the ‘direct* problem of spherical geometry, which is: Given 
the latitude, and longitude. A, of a point Fi(^i,Ai), the distance i and the forward 
asimnth ai 3 to a point i^(^, A 3 ), determine A 3 and the backward asimuth 031 . 

Solve for ^ and A 3 . Relabel the angles and legs so that A = 013 , B — X\- A 3 , 
C = 360® - «i 3 , a = 90® - 6 = d and c = 90® - ^ 1 . Substituting these into Equs. ( 1 ) 

and simplifying the differences of angles, we obtain 

sin ^3 = cosdsin^i +sindco 8 ^i cos 0 x 3 , 
cos^cos(A 3 — Ax) = cosdcos^x -sindsin^xco 8 ax 3 , 
cos^sin(A 3 - Ax) = sindsinax 3 . 

The first equation can be used to solve for Using the principle angle of the arcsine, ^ 
lies in the range —90® to 90®, which is correct for a latitude. Since ^ lies between —90® 
and +90®, the sign of cos ^ is known, thus the system has only one solution. Further, and 
again most important, the sign of cos ^ is positive for all ^3 between —90® and +90® and 
so there are no special cases. 

The first equation is used to compute 

^ = arc 8 in(cosdsin^x +flindcos^x cosax 3 ). ( 6 ) 

The second and third equations determine cos(A 3 - Ax) and sin(A 3 - Ax). As before, with 
both the sine and the cosine known, the quadrant is uniquely determined using the qatn 
function. Since cos ^ is always positive, it can be eliminated in both equations. Then 

A 3 = Ax +qatn(sindsinax 3 ,co 8 dco 8 ^x — sindsin^x cos 0 x 3 ). (7) 

The back asimuth can be found by relabeling the angles and legs so that A = 03 x, 
B = 360® — flt 3 x, <7 = Ax — A3, 0 = 90® — 6 =s 90® - ^x and c = d. Substituting these 

into Equs. ( 1 ) and simplifying the differences of angles, we obtain 

0 f 3 x =qatn(— coa^x 8 hi<»i 3 , 8 iu^x 8 ind — cos^xcoodcosaxs). ( 8 ) 
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This second fonnula for the back asinmth is need in the spirit that the unknown parameters 
shonld be determined as a fonction of the known parameters rather than determining 
the second unknovm parameter as a fonction of the previously evaluated first unknown 
parameter. That is, any one of the unknowns can be found as a function of the knowns 
without having to first find one of the other unknowns. 

Hhutrmtlon S. Although derived from different principles, the equations of coor- 
dinate transformation obey the rule of Chauvenet. ibr example, the equations which 
transform local hour angle and declination to altitude and aiimnth are 

cosasin A = — cos^sink, 

cosdcosA = sin^cos^ — cos^sin^cosk, (9) 

sina = sin^sin^ + cos^cos^coek, 

where a = altitude, A = asimuth, S = declination, h— local hour angle and ^ = latitude 
[Bef. 4, pg. 26|. 

The third equation is used to solve for altitude 

a = arcsin(sin58in^ + co8^co6^co8k). (10) 

The altitude lies between -90^ and +90*’, and so is correctly determined by the arcsin 
function. Since cos a is always positive it can be eliminated between the first two equations 
when using the qatn function. The local hour angle, which lies in the range 0^ to 360^, is 
determined by 

A = qatn(— co8^8ink,sin5co8^ — cos^sin^cosh). (11) 

C. Simplification of Standard Works. Using the general triangle and the qatn func- 
tion, the simplifications listed below can be made in Volume II of Bowditch [Ref. 1]. The 
various sections are preceded by the section symbol § and equation numbers from Bowditch 
within those sections are prefixed with the letter B. Unprefixed equation numbers are those 
contained in this document. 

§700. Solving for altitude. — Using the conventions that northern latitudes and 
declinations are positive and that southern latitudes and declinations are negative, the 
altitude of a star can be computed directly from Eqn. B(2a) or Equ. B(2b). Equ. B(2a) is 
equivalent to Equ. (10). All of the special cases are thus eliminated. 

§707. Solving for asimuth.— Instead of using Equ. B(4b) or Equ. B(5b), use 
Equ. (11). See §708 below. 
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§708. Cnme ummth. — Bowditch uses meridian angles which can be labeled either 
east or west. Since altitnde/aiimnth is a left-handed coordinate system, the convention is 
that west meridian angles are positive and east meridian angles are negative. 

Example i.— The latitude of the observer is 30^ 26/ON; the declination of the celestial 
body is 22*06/2N; the local hour angle is 38"64/7W. Using Equ. (11), 

A = qatn(-cos22*06/28in39*54/7,sin22®06/2co630*25/0 

- cos22*06/2sin30®25/0cos39®54/7) 

= qatn|-(0.92661)(0.64161), (0.37628) (0.86237) 

- (0.92651)(0.50628)(0.76703)] 

= qatn(-0.59446, 0.32449 - 0.36980) 

= qatn(-0.69446, -0.03631) 

= -93.®39913 
= 266.®60087 
il = 266®36/l. 

Example 2, — ^The latitude of the observer is 30^ 25/OS; the declination of the celestial 
body is 22*06/2N; the local hour angle is 39^64/7E. Using Equ. (11), 

A = qatn(- cos22®06/2 sin(-39* 54/7), sin 22®06/2cos(-30*25/0) 

- cos22®06/28in(-30®25/0) cos(-39®54/7)) 

= qatn(-(0.92661)(-0.64161), (0.37628) (0.86237) 

- (0.92651)(-0,50628)(0.76703)] 

= qatn|0.69446, 0.32449 - (-0.36980)) 

= qatn(0.69446, 0.68429) 

= 40.®98141 
A = 40^58/9. 

D. Sphwoid Eartli and Great Circle Fommlas. 

1. Spheroid Earth Fomralas. The spheroid earth algorithms are those of Thomas 
[Ref. 5|. Some changes by Shudde [Ref. 6] make quadrant determination automatic by 
use of the qatn function. Eiast longitudes and south latitudes are negative. The earth’s 
equatorial radius is denoted by a« and the earth’s flattening factor is denoted by /. 
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Sk lnTU06 SolntioB* Lipnt vaiiabltt 3n the latHiide and longitude Ai of Pij 
the latitude ^ and longitude h of Pi, h ^<1 Output variables are distance, (Sja,), 
foTward asimnth hrom P\ to Pi, an, and bade asimnth from Pi to Pi, an. All an g ular 
input and output is in radians. Compute: 

fi = tan"‘ [(1 - /) taadil, for • = 1, 2, 

AA = A, - Ai, tfm = (»i + »a)/2, = (fi - fi)/2, 

H = cos* Mm - siu* im, i = shi* Mm + B 8in*(AA/2), 
i = 2sin-‘(i‘/*), tr = 2rin* cos* Ad„/(1 - L), 

7 = 2sin*A«„coe*tf„/i, JT-tf + V, Y = U-V, 

T = i!wii, Z1 = 4T*, S = 2cosd, A = DE, 

C = T~(A-E)I2, ni = X{A + CX\, 

B = 2D,m = Y{B + EY),m = DXY, (12) 

iid = f{TX- Y)/A, Sid = /*(n, -m + nj)/64, 

5/a, = (r-«ii + «jd)8ind, M = Z 2 T - [20T ~ A)X - (B + *)Y, 

F = 2Y-E{4-X),G = /T/2 + /*Af/64, Q = -(f Gtan AA)/4, 

AAi, = (AA + (J)/2, 

<1 = qatn(— 8inAdmCosAA{„,coedm8inAAj„), 
t, = qatn(cos Af„ cos 4A;„, sin sin AA(„), 

Ofia = ti + and ofQi — ti — ia. 

b. Direct Sohition. Input variables are the latitude <f>\ and the longitude Ai of 
Pi, the forward arimuth ocia from Pi to Pa, and the distance (5/a«) from Pi to Pa. The 
output variables are the latitude ^ and longitude of P2z and the backward azimuth 
a^i from P3 to Pi. All angular input and output is in radians. Compute: 

=tan"^[(l-/)tan^i], 

M = co8^ifflnai3, ci = /Af, ca = /(I- Af^)/4, 

D = {l — ca)(l — Ca ” ciAf), P = ca(l + CiM/2)/D, 

N = cos^i cosaia, (Ti = qatn(Af,8in9i), 
d = {S/a^)/Df tt = 2 ((Ti - d), FK = 1 - 2Pcosu, 

V = cos(tt + d), X = c? sin dcosd(2K^ - 1), 

r = 2P7W’sind, A(T = d + X-r, (13) 



8 



K = [AP + {NcosAa - am^i an 
tan ^3 = (sin^iC08A(T4*iV^8mA(r)/ir, 
h = tan*M(tan^3)/(l - /)!, 

Ati = qatn(flin AiTflin ai3, cos cos A<r - sin sin A<r cos 013), 

H = ci(l - C3 )Ao’ - C1C3 sin A<tcos( 2(71 - A(t), 

A3 = Ai + Atf — Hf 

ctii = qatn(-Af, ~{N cos A(7 — sin sin A(t)]. 

2. Spherical Earth Fomiiilaa. The inverse solution [Equs. ( 3 ), ( 4 ) and ( 5 )) and direct 
solution [Equs. (6), ( 7 ) and (8)] formulas were developed in Section B of this report. They 
are summariBed here for convenience. 

a. InverM Solution* Input variables are the latitude and lon^tude Ai of Pi, the 
latitude ^ and longitude A3 of P3, and A3. Output variables are distance, d, forward 
asimuth horn Pi to P3, 013, and back asimuth from P3 to Pi, CE31. All angular mput and 
output is in radians. Compute; 

d = arccos|sin ^3 sin ^1 + cos ^3 cos ^1 cos(A3 - Ai )), 
oria =qatn(co8^sin(A3 - Ai),sin^cos^i -cos^sin^iCOs(A3 - Ai)), and ( 14 ) 
031 =qatn(-cos^isin(Ai - A3),sin^i cos^ -cos^isin^cos(Ai -A3)). 

b. Direct Solution. Input variables are the latitude ^1 and the longitude Ai of Pi, 
the forward asimuth 013 frx)m Pi to P3, and the distance d frnom Pi to P3. The output 
variables are the latitude ^ and longitude A3 of P3, and the backward asimuth 031 from 
P3 to Pi. All angular input and output is in radians. Compute: 

^ = arcsin(cosdsin<^i +sindco84ii cosoi3), 

A3 = Ai +qatn(8indsinoti3,cos(icos^i — sinifsin^icosais), and ( 15 ) 
< 3(31 = qatn(— cos ^1 sin 013, sin ^1 sin d — cos ^1 cos d cos 013). 

c. Given Longitude) Find Latitude* Input variables are the latitude ^1 and 
longitude Ai of Pi, the forward asimuth 013 from Pi to P3, and the longitude A of some 
point P on the great circle joining Pi and P3. The output variable is the latitude of 
P. All angular input and output is in radians. Using Equs. (2), label the angles and legs 
so that A — 360 ® — 031, B = ori3, (7 = A — Ai, a = 90 ® — ^1, 6 = 90 ® - ^ and c = d. 
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Sttbetituting these into Eqna. (2) and simplifying the differences of angles, we obtain 

costt 2 i = — costtiaco^A — Ai) + 8maiasin(A — Ai)sin^i, 
sinoai sin^ = cosaiasin(A — Ai) + sin oia cos(A — Ai)sin^i, and 
sinaai cos^ ~ ^aiacos^i. 



Since the only quantity of interest is and since latitude lies in the range of —90® to +90® , 
it suffices to divide the second equation by the third equation to obtain an expression for 
tan We obt^ 



^ = tan”^ 



' co8aiasin(A - Ai) +8inaiacos(A - Ai)sin^i ' 
sin aia cos 



(16) 



Fbr almost all applications, Equ. (16) will suffice. It is, however, possible for the denom- 
inator of Equ. (16) to become sero in certm circumstances. We now face a situation in 
which special cases must be examined. 

One way to handle the dilemma is to first evaluate the denominator and the numerator. 
K the value of the denominator is sero, then ^ 90® if the numerator is positive or 

^ = -90® if the numerator is negative. If the d^ominator is not sero, then use Equ. (16). 
A second way to handle the dilemma is to compute 



^ = qatn(cosoia sin(A - Ai) + sin aia cos(A - Ai) sin ^i, sin aia cos ^i], (17) 

but this leads to another problem because ^ may now be in the range of - 180® to +180® . To 
correct the output to a proper latitude, subtract 180® if ^ > 90® or add 180® if ^ < —90®. 

The question of which method to use is a matter of personal choice. Generally, a 
simple one line logic function can handle all of the special cases. 

d. Find the Vertex. The vertex of a great circle route is the most northerly or/and 
the most southerly point on the great circle. Bowditch [Ref. 1, §1016, Equs. B25-B30] gives 
formulas for first, computing the latitude of the vertex, and second, computing the longi- 
tude of the vertex as a function of the latitude of the vertex. The Bowditch latitude formula, 
from a straight forward ^plication of the Law of Sines, is = cos~^ (sin oia cos ), where 
is the latitude of the vertex and and oia have the same meaning as in the previous 
section. The principle angle of the arccosine lies between 0® and 180® whereas latitudes 
lie between —90® and +90®. Bowditch illustrates only cases in which is between 0® and 
+90® and, unfortunately, makes no mention of how to handle cases when the arccosine is 
greater than 90®. 
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There are several alternate methods of developing equations for the vertices. One 
method is presented here and an alternate method is presented at the end of §e. Input 
variables are the latitude and longitude X\ of Pi, and the forward asimuth ais from Pi 
to P 3 . Let the latitude and longitude of some arbitrary point P on the great circle from 
Pi to Ps be denoted by (j> and A, respectively. All angular input and output is in radians. 
Label the angles snd legs so that A ~ A — Ai, P ~ aia, C ~ 360^ — 0 ( 31 , a = — 90®-^ 

and e = 90” - ^ 1 . Substituting these into Equs. ( 1 ) and simplifying the differences of 
angles, we obtain 



cosd = sin^sin^i +cos^cos^i cos(A - Ai), 
smdco 8 ai 3 = sin^co 8 ^i -cos^sin^icos(A-Ai), and 
sindsinai 3 = cos^sin(A- Ai). 



Divide the third equation by the second equation to obtain 

cos^ 8 in(A- Ai) 



tanai3 = -r- 



sin^cos^i — cos^sin^i cos(A — Ai)' 

Then rearrange the equation in the form 

tanai3[sin^cos^i ~co8^sin^iCos(A - Ai)] =co8^sin(A - Ai). (18) 

Next, we wish to find what value of ^ is an extremum of A. To do this we must first 
compute d<^/dX in Equ. (18): 



tanai3 



cos^cos^i^+sin^sin^i cos(A - ^ 1 )^ + cos^sin^i sin(A - Ai) 



SB — sin^ 8 in(A — ■^ 1 )^ + co 8 ^cos(A - Ai). 



Then we must set 



Thus, 



dX 



4 - 4 ^ - 



tan ai 3 cos sin <^i sin(A« - Ai) = cos ^ cos(A« - Ai), 
which rearranges to 



or 



tan(A* - Ai) = ^ 
A, = Ai + tan“^ 



1 



sin<^i tanai 3 
1 



sin^i tanai 3 



(19) 
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Once A« is known, can be found using one of the methods in §c such as Equ. ( 17 ). Since 
|A« - All < 90 **, it follows that ^ must be in the same hemisphere as Thus, if ^i is in 
the northern (southern) hemisphere, then must be the northern (southern) vertex. 

e. Given Ladtnde^ Pind Longitude. Input variables are the latitude <j>i and 
longitude Ai of Pi, the forward asimnth ai^ hrom Pi to P] and some latitude It is 
convenient but not necessary to know the latitude of the vertex of the great circle hx>m 
Pi to P). We wish to find the longitude(8) on the great circle from Pi to P3 which have 
latitude There are several possible outcomes: (1) if |^| > |^«|, there is no solution; (2) 

if ^ there is one solution and that is A = A«; and ( 3 ) if |^| < |^«|, there are two 

longitudes which satisfy the given conditions. 

Rewrite Equ. ( 16 ) in the form 

cos aia cos ^ sin(A - Ai ) + sin aia sin ^1 cos ^ cos( A - Ai ) = sin aia cos ^1 sin 



or 

5 sin(A - Ai) + Ccos(A - Ai) = if, (20) 

where S = cosancos^, C = sinaiasin^i cos^, and K = sinaia cos^i sin^. Define p >0 
and tf as the solutions to the equations 

C = pcosrj and 5 = />sinf;. (21) 

Solving, we find that 

p = +>/52 + C^ and 
ly = qatn( 5 , C) 

Substituting Equs. (21) into Equ. (20) and rearranging, we find 

cos(A - (Ai + »;)] = Kip, 



So that 



A-(Ai + iy) = ±cos-i {Kjp), 



or 

A = (Ai + fj) ± cos"^ [KIp) . (22) 

Recall that a knowledge of is not known to determine the number of solutions. (1) if 
\Kjp\ > 1, there are no solutions; (2) if \Kjp\ = 1, there is one solution; and (3) if 
\K/p\ < 1, there are two solutions. 
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In the previous section it was mentioned that there was yet another way to detennine 
the position of the vertex. To do this, for fixed define /(A) in terms of £!qii. (20) as 

/(A) = 5sin(A - Ai) + Ccos(A - Ai) - K, 

Then /'(A), the derivative of /(A) with respect to A is 

f(X) = 5cos(A - Ai) - Csm{X - Ai) 



Define A* so that /'(A*) = 0. Thus 





5cos(A® - Ai) - Csin(A® - Ai) = 0, 


or 


sin(A*-Ai) S 
cos(A®-Ai) ■ C" 


or 


A* - Ai = qatn(5, C), 


so that 


A* = Ai + qatn(5, C) = Ai + ty. 



Hence A* = Ai + is an extremnm of Equ. (20), or A« = A*. Note that 

/«(A*) = -5sm(A* - Ai) - Ccos(A* - Ai) 

= — fifsinry — Ccosiy 
CP 



P P 
S^+CP 
P 

= -p < 0. 



Thns Ay = A'* is the northern vertex. 

S. Rhnmb Line (Mercator) Pomralas. 

1. Spheroid Sarth Fommlas. Bowditch [Ref. 1, Explanation of Ihble 5, Meridional 
Parts] gives the formula for the computation of meridional parts as 

M = a« Intan As® + 






= Uelntan ( 45® + 



■ 






I «sm^+ ySl 



sin^^ + T + 
5 



■) 
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where M is the number of meridional parts between the equator and the given latitude 
Oe is the equatorial radius of the earth, and « = y/[2f - is the eccentricity of the 
earth, where / = 1/298.26 is the WGS 1972 earth flattening factor. 

Using the id^tity |R^. 7, #769] 







I*" <11 



it is possible to rewrite the formula for M as 



M= ..liitaii (4$« + 1) - 

la tan ^45'> + - la 



= «« 



= a«ln 



/l + esin 
^1 — esin^ j 



If we denote the course (forward asimuth) by a^, then [Bowditch, VoL n, §1001 and 
§1013] AA = Aj - Ai = mtanoia where m = M(^) - M{^\)y so that 



tan ai <2 = 



tan 



In 



(A2-Ai)/a« 

_ taa (4S» + 



(23) 



/ I + esin<^ / 1 + eain^i 

esin^y Y-esini^iy 



This equation is also given by Shufeld and Newcomer [Ref. 8 , pp 81-84]. It is the spheroid 
earth rhumb line equation. 

a. Inverse Solution. Input variables are the latitude and longitude Ai of Pi , the 
latitude <j>^ and longitude A^ and P^, ^ and A). Output variables are distance, d, forward 
asimuth ai^ from Pi to P^, and back asimuth a^i from P 3 to Pi. All angular input and 
output is in radians. In Equ. (23), express the angles in radians, then to obtain 013 in the 
proper quadrant, rewrite Equ. (23) in the form 
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The distance is 



<»< i 



if cosaia ^ 0 , 



d=l "'cosaia' ' (25) 

a«|Aa — Ai I cos , otherwise. 

In the WGS 1972 earth model, a« = 3443.9174 nautical miles. Also, aai = aia + r. 

b* Direct Solution, bipnt variables are the latitude <f>i and the longitude Ai of Pi, 
the forward asimuth aia from Pi to and the distance d from Pi to Pa. The output 
variables are the latitude ^ and longitude Aa of Pa, and the backward asimuth aai from 
Pa to Pi. AH angular input and output is in radians. In the direct rhumb line solution for 
the spheroid earth, if cosaia = 0 , then 



Aa = Ai + 



cos^i 






(26o) 



otherwise, if cos aia 



Aa = Ai + tan aia 



/ I + esm^a Y 
_ 1 1 — esin^ j 
^ = ^i + (d/a*)cosaia. 




/ I -t-g^^i 
\l - esin^i j 



, and 



(266) 



Also, aai ~ <*13 + **• 



3. Spherical Earth Fommlas. The spherical earth fonnulas are obtained by setting 
e = 0 in Equations (24) and (26b), and replacing Oe by the standard approximation of 60 
nautical miles po* degree of arc on the earth’s surface. 

a. Inverse Solution, hiput variables are the latitude <ffi and longitude Ai of Pi, the 
latitude ^3 and longitude A^ and P 2 . <l >2 and A^. Output variables are distance, d, forward 
asimuth from Pi to P 3 , ai^, and back asimuth from Pg to Pi, a^i. All angular input and 
output is in radians. 

The spherical earth rhumb line course is obtained by setting e = 0 in Equ. (24). The 



result is 

aia = qatn j^Aa - Ai, lntan - In tan 

For the sph^cal earth model, the distance d in nautical miles is 

_ J \r J cosaia 
’ 60 



d = 






•) - 



Ai|cos<^i, otherwise. 



(27) 



(28) 
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Also, otai = aria + k. 

b. Direct Solntion. bipnt variables are the latitude and the longitude Ai of Pi, 
the forward aiimuth aia from P\ to Pa, and the distance d bx>m Pi to Pa. The output 
variables are the latitude ^ and longitude Aa of Pa, and the backward asimuth aai from 
Pa to Pi. All angular input and output is in radians. 

If cosoria — 0, then 

^ (lio) (m) ^ (29a) 

<h = hi 

otherwise, if cosaia ^ 0, 



^2 = + 

Also, fltai = Qfia + s’. 



[■"‘“(f + y) 


-•"‘“(j+t)] 


(^) 





, and 



(296) 
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m. NAVALGOR: Navigation Algorithm Program. 

A. Introduction. NAVALGOR implements three sets of computational procedures: (1) 
Direct Solution Algorithms, ( 2 ) hiverse Solution Algorithms, and (3) Rhumb Line Ap- 
proximations to Great Circle Routes. 

The *Direct Solution Algorithms* compute the latitude and longitude of a position P 3 
and the backward asimuth from P 3 to Pi given the latitude and longitude of a position Pi, 
the forward asimuth from Pi to P3 and the distance from Pi to P3. The “Inverse Solution 
Algorithms* compute the distance from position Pi to position P 3 , the forward asimuth 
from Pi to P 3 , and the backward asimuth from P 3 to Pi given the latitude and longitude 
of positions Pi and P3. For comparision, the direct and inverse computations are made 
using four procedures: ( 1 ) the spheroid earth model, ( 2 ) the spherical earth model, and 
rhumb line approximations using both (3) the spheroid earth and (4) the spherical earth 
models. 

The spheroid earth model should be considered the standard of comparison for all 
other models. Computations performed using an IBM soss and double precision FORTRAN 
to compare the inverse solution algorithm [Equs. (12)] to the AGIO long lines (50 to 6000 
miles) (Ref. 9] demonstrated that the maximum average error of the inverse solution algo- 
rithm was -0.005 meters with a standard deviation of 0.014 meters (unpublished result). 

The spherical earth model is included because of its popularity and simplicity. For 
many procedures the spherical earth model is adequate, and it certainly is more compact 
for use in a small computer than the spheroid earth model 

The “Rhumb Line Approximations to Great Circle Routes* procedure may be used to 
find piecewise constant course (rhumb line) approximations to great circle routes for any 
given increment in longitude. In addition, if the vertex of a course is, for example, too far 
to the north, a limiting latitude may be input to restrict the rhumb line approximation to 
go no further north than that limiting latitude. The spheroid earth rhumb line equations 
are used in this section of the program. 

In the sample problems and in the program listing, the east minus and south minus 
convention is used. 
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B. Sample Problema — ^NAVALGOR# The output shown in these sample problems 
was computed using Microsoft BASIGA/D (double precision) on an IBM PC. 

Maater Mcira. The master menu for NAVALGOR is shown below. 

NAVIGATION ALGORITHM DEMO 

1) DIRECT SOLUTION 

2) INVERSE SOLUTION 

3) RHUMB LINE APPROXIMATIONS 

TO GREAT CIRCLE ROUTES 

4) QUIT 

Option (1) selects the "Direct Solution Algorithm*, option (2) selects the "Inverse Solution 
Algorithm*, option (3) selects the "Rhumb Line Approximations to Great Circle Routes* 
algorithm, and option (4) returns the user to the operating system. 

Problem 1. Suppose you are at San fVancisco (latitude 37** 47' north and longitude 
122"25' west), that your initial course is 260" and that you travel a distance of 4000 n. mi. 
What is your final position? Select Option 1 from the master menu. 

DIRECT SOLUTION 

1st LATITUDE DD.MMSS (-S) 7 37.47 
1st LONGITUDE DDD.MMSS (-£) 7 122.25 
INITIAL COURSE DDD.MMSS 7 260 
DISTANCE (n. nl.) 7 4000 

SPHEROID EARTH DIRECT SOLUTION 
2nd UTITUDE 6"38.7' 

2nd LONGITUDE -172"08.8' 

BACK AZIMUTH 51"40.7' 

SPHERICAL EARTH DIRECT SOLUTION 
2nd UTITUDE 6"41.9' 

2nd LONGITUDE -172"00.7' 

BACK AZIMUTH 51"36.9' 

RHUMB LINE SOLUTIONS 

SPHERICAL 26"12.4' -159"66.0' 

SPHEROID 26® 12. 4' -160® 18. 4' 

PRESS ANT KEY TO CONTINUE 
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Problem 2. Suppose you are at San fVancisco (latitude 37” 47* north and longitude 
122”25' west) and that your destination is Sydney, Australia (latitude 33”5P south and 
longitude 151”13^ east). How far do you travel, what is your initial course, and what is the 
backward asimuth from Sydney to San Fransisco? Select Option 2 from the master menu. 

INVERSE SOLUTION 

let LATITUDE DD.MMSS (-S) 7 37.47 
1st LONGITUDE DDD.MNS8 (-E) 7 122.25 
2nd UTITUDE DD.MMSS (-S) 7 -33.61 
2nd LONGITUDE DDD.MMSS (-E) 7 -151.13 

SPHEROID EARTH INVERSE SOLUTION 
DISTANCE 6443.52 n.isi. 

FORWARD COURSE 240”29.3' 

BACK COURSE 65”65.7' 

SPHERICAL INVERSE SOLUTION 

DISTANCE 6446.3 n.al. 

FORWARD COURSE 240”18.9' 

BACK COURSE 65”45.9' 

RHUMB LINE SOLUTIONS: SPHERE SPHEROID 

DISTANCE 6464.49 6485.71 n.al. 

COURSE 228”19.7' 228”29.7^ 

PRESS ANY KET TO CONTINUE 

Problem S. Suppose you are at latitude 37” north, longitude 76” west and that your 
destination is latitude 45” north, longitude 1” west. Compute the initial great circle course 
and distance, the latitude and longitude of the vertex, and a rhumb line approximation to 
the course traveling at most 7” degrees of longitude on each leg of the course. In addition, 
no portion of the course is to be more northerly than 47”. Select Option 3 from the master 
menu. 
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iBt Screen: 



RHUMB LINE APFROXINATIONS 

INITIAL LATITUDE DD.MMSS (-S) 7 37 
INITIAL LONGITUDE DDD.MNSS (-E) 7 76 
FINAL UHTUDE DD.MMSS (-S) 7 46 

FINAL LONGITUDE DDD.MNSS (-E) 7 1 

GREAT CIRCLE SOLUTION 

INITIAL COURSE 56^21.2' 

DISTANCE 3307.8 n.oi. 

RHUMB LINE (MERCATOR) SOLUTION 
COURSE 81<^58.2' 

DISTANCE 3436.9 n.ol. 

VERTEX: LAHTUDE 48®19.8* 

LONGITUDE 28®07.3* 

PRESS ANT KEY TO CONTINUE 

In this first screen, the great circle solution and a single riiumb line solution are computed 
and displayed. In addition, the location of the nearest vertex is displayed. The vertex 
may or may not lie between the initial and final positions. If the vertex is not between 
the initial and final positions then any limiting latitude which may be input on the next 

screen prompt will be ignored. If the vertex is between the initial and final positions, a 

limiting latitude will be ignored unless it is lies between the latitude of the vertex and the 
latitude of the position which is closest to the latitude of the vertex. 

2nd Screen: 

RHUMB LINE APPROXIMATIONS 

INPUT THE MAHMUN NUMBER OF DEGREES 
OF LONGITUDE ON EACH RHUMB LINE LEG 7 7 

LIMIHNG LATITUDE DD.MMSS (-8) 

(0 TO OMIT) 7 47 
LIMITING LONGITUDES: 

10® 46. S' 46®28.8' 

PRESS ANY KEY TO CONTINUE 

In the second screen a prompt appeared for the maximum longitude between each rhumb 
line leg — a value of 7® was input. Also input was a limiting latitude of 47®, which lies 
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between the vertex latitude of 48*^19.8' and 45^, the latitude of the final position. Limiting 
longitudes of 10*45.8' and 45*28.8' are computed and displayed. The latitude of a great 
circle route between these limiting longitudes will be more northerly than 47* , consequently 
any thumb line approximation will be due east or due west at a latitude of 47* between 
these longitudes. 

3rd Screen: 



RHUMB LINE APPROXIMATION TO GREAT CIRCLE COURSE 







GREAT 


GREAT 


RHUMB 


RHUMB 






CIRCLE 


CIRCLE 


LINE 


LINE 


UT 


LONG 


COURSE 


DIST 


COURSE 


DIST 


37*00.0' 


76*00.0' 


56 


0 


58 


0 


39*27.9' 


71*00.0' 


69 


278 


62 


279 


42*18.8' 


64*00.0' 


64 


639 


66 


641 


44*32.0' 


57*00.0' 


69 


971 


71 


974 


46*11.7' 


60*00.0' 


74 


1283 


76 


1287 


47*00.0' 


46*28.8' 


77 


1476 


90 


1480 


47*00.0' 


10*46.8' 


103 


2884 


106 


2900 


46*56.3' 


6*00.0' 


107 


3130 


108 


3148 


45*00.0' 


1*00.0' 


110 


3308 


108 


3326 



1) HEN APPROXIMATION OR 2) NEW PROBLEM? 

Screen 3 displays the rhumb line approximation for a TnaximuTn of 7* of longitude between 
course changes. The latitude and longitude of the initial and final positions and of the 
positions at which course changes are made are given in columns 1 and 2, respectively. 
Column 3 shows the great circle heading at each position, while Column 4 shows the 
cumulative great circle distance from the initial position. Similarly, Column 5 shows the 
rhumb line course to be followed between each pair of positions and Column 6 shows 
the cumulative rhumb line distance from the initial position. Note that at 47* north and 
45*28.8' west (the first limiting latitude position) the rhumb line course changes to due east 
until 47* north and i0*45.8' west in reached. Also note that the rhumb line approximation. 
Including the Metour' at the limiting latitude, is only 18 n. mi. longer than the great circle 
route. The third screen ends with two options: Option 1 returns to Screen 2 and prompts 
for new inputs while Option 2 returns to the master menu. 
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C* Program Lktlng — ^NAVALGOR 



10 SOI "NATIGiTION ALGOEITBNS* ft.H. SHODOE, 03*12-86. BET. 00-14-86 1600 
13 BEK "NAVALGOB/BAS* 

30 OEFDBL A-Y 

40 P4»ATH(1) :P>P4+P4:P>P2+P2:TP-PI+PI:BD-PI/180;EP«lE-33 
60 FL-l/208.26:AE-6378138/1862:SC-SqR(FL*(2-FL)) 

60 DEFFI«a)-X-MO«INT(X/llO):BEK X MOD NO FUNCTIQ}!. 

70 DEFFHL(I)-X-TP*I!IT((I+PI)/TP):BEK LOMCITUDE ADJUST (-PI.PI) 

80 0EFFNR(X)-INT(X«N0+.6)/H0:BEK BOUNDING FUNCTION. 

00 0EFFNG(X)-X<»PI*8GN(X)*(AfiS(X)>P2):BEM LATITUDE ADJUST (-PI/2.PI/2) 

100 OEFFNC(Z)-ATN(sqB(l-X«X)/(X-EP«(X^)«)))-PI*(X<0«):BEN ABCCOS 
110 D£FFN8(X)»ATN(X/(SqB(l-X*X)-EP*(AB8(X)»l))):BQI ABCSIN 
116 GOTO 2000 

120 BE3f QATN (-PI.PI) FUNCTION: 

130 A»AIN(Y/(X-EP*(X= 0 #)))-PI*(X< 0 #)*(SGN(T)-(Y= 0 #)) :aETORN 
140 BEN QATN (O.TWOPI) FUNCTION: 

160 A»AIN(Y/(X-EP*(X=0#)))-PI*(X<0#)+TP*(X>»0#)*(T<0i) :BETUBN 
170 BOi 

200 BDI DIBECT SOLUTION. SPHEBOID EABIH. ALL ANGLES NUST BE IN RADIANS. 

210 BOI INPUT: UTITUDE Gl. LONGITUDE LI. FOBVABD AZIMUTH A1 AND 
220 BOI DISTANCE OD-(S/AE) TO A POINT P2. NOTE: DD HAS RADIAN UNITS. 

230 BOI OUTPUT: UTITUDE G2. LONGITUDE L2 AND BACXWABO AZIMUTH A2. 

240 8d-8IN(Al) :C0-C0S(A1) :TA-ATN((1-FL)«TAN(G1)) 

260 S8«8IN(TA):C8-C08(TA):N— S0*C8:C1-FL*M 
260 C2-FL*(1-M*M)/4:D-(1-C2)*(1-C2-C1*M) 

270 P-C2«(1<^.6«C1*N)/D:N-C8«C0 

280 y^N:X-S8:G08UB 13O:8G-A:SO-0D/O:U-2*(SG-SO) 

200 N^1-2*P*C08(U) :V=C08(U+8D) :X=C2*C2*8IN(8D)*C08(8D*(2*V*V-1)) 

300 Y^«P*V*W«8IN(SD) :Dq-8D^X-Y:SL-8IN(D0) :CL-COS(Dq) 

310 X-8qB(M*M+(N*CL-88*8L)'‘2) :G2-ATN((88*CL+N*8L)/X/(1-FL)) 

320 Y— SL*S0:X-C8«CL-S8*SL*CS:G0SUB 130:DJ-A 
330 H-Cl*(l-C2)*D9-Cl*C2«8L«C0S(SG^86-0q) 

340 L2=FHLai+DJ-H):Y»M:X=-(N*CL-88*8L):G08UB 150:A2=A 
360 BETUBN 
360 REM 

400 REM INVERSE SOLUTION, SPHEBOID EARTH. ALL ANGLES MUST BE IN RADIANS. 
410 REM I NPUT: UTITUDES Gl A G2. AND LONGITUDES LI k L2. 

420 BEM OUTPUT: DISTANCE DD»(S/AE), FOBWARD AZIMUTH Al, AND BACKWARD 
430 REM AZIMUTH A2. NOTE: DD HAS RADIAN UNITS. 

440 TA-ATN((1-FL)*TAN(G1)) :TB-ATN((1-FL)*TAN(G2)) 

460 DM-L2-Ll:DT»(TB-TA)/2:TM-(TA+TB)/2 

460 SH-SIN(DT) :CH-COS(DT) :SM-SIN(TM) :CM-COS(TM) 

470 H-CH*CH-SN*8M:L«8H*8Hi>H*(SIN(DM/2))''2 

480 SD^*FN8 (SqR(L) ) : U-2*SM* SM*CH*CH/ ( ( 1 -L) - (EP* (L-1) ) ) 

400 V«2*8H*8H*CM*CM/L:X»U+V:Y-U-V:>8D/8IN(8D) 

500 D-4*T*T:E-2*C08(8D) :A-D*E:OT-(A-E)/2:Nl-X*(A+C*X) 

610 B-D+D:N2-Y*(B+E*Y) :N>D*X*Y:D4»PL*FL*(Nl-N2+N3)/64 
620 D3»FL*(T*X-Y)/4 :DD-(T-D3+D4)*8IN(8D) 

630 M»32*T- (20*T-A) *X- (B+4) *Y : F=Y+Y-E*(4-X) 

640 G-FL*(T/2+FL*M/64) :q—(F*G*TAN(DM))/4 
660 DW-(DM-^q) /2 : SW-SIN(DW) : CW-COS(DW) 
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660 

670 

680 

600 

600 

610 

620 

630 

640 

660 

660 

670 

680 

600 

700 

710 

720 

730 

740 

750 

760 

770 

780 

700 

800 

810 

820 

830 

840 

850 

860 

870 

880 

885 

800 

000 

010 

1000 

1010 

1020 

1030 

1040 

1050 

1060 

1070 

1080 

1000 

1100 

1110 

1120 

1200 

1210 

1220 



T>8H*CV:X-GN«8W:6a8UB 130:T1«A 
Y«'CB*GV:X«8N*8W:6a8DB 130:T2«A 
MO-TP: A1-PI0I(T1+T2) :A>FID1(T1-T2) :RE7Umf 
BEN 

BEM DIBECT SOLDnON. SPHEBIGAL EABTH. ALL ANGLES MUST BE IN BADIAN8. 
BEN INPUT: UnTDDE Gi. LONGITUDE Li, FOBMABD AZIMUTH A1 AND 
BEN DISTANCE DD TO A POINT P2. NOTE: DD IS IN BADIANS. 

BEN OUTPUT: UTITUDE G2. LONGITUDE L2 AND BAGEWABD AZIMUTH A2. 
SO-SIN(Al) :C0=^:08(A1) :S3>8IN(Gi) :C3«C08(Gi) 

SD-8IN(DD):CD-C08(DD) 

T>8D*8O:X-C3*C0-83*CO*8D:G08UB 130:L2»FNL(L1-A) 
T>-S0«C3:X«8D*83-CD*C0*C3:G0SUB 160:A2>A 
G2-FN8(S3*CD+C3*8D*C0) : BETUBN 
BEN 

BEN INVEBSE SOLUTION. SPHEBIGAL EABTH. ALL ANGLES MUST BE IN BADIANS. 
REM INPUT: UTITODES Gi k G2. AND LONGITUDES LI k L2. 

BEN OUTPUT: DISTANCE DD TO A POINT P2. (NOTE: 0 <> DD <» PI BADIANS). 
BEN POBHABD AZIMUTH Al. AND BACXWABD AZIMUTH A2. 

S3-8IN(G1) :C3-C08(Gi) :84-8IN(G2) :C4-C08(G2) 

DL»L1-L2:SU>8IN(DL) :CU=C08(DL) 

DD-PNC(S3«84-K:3*C4«CU) 

>8U*C4:X-C3*84-83«C4*CU:G0SUB 150:A1-A 
T»8U«C3:X«C4«83-84*C3*C0:G08UB 160 :A2«A: BETUBN 
BEN 

BEN INVEBSE SOLUTION. BHONB LINE FOB SPHEBIGAL k SPHEROIDAL EABTH. 

BEN INPUT: UTITUDES Gi k G2. AND LONGITUDES LI k L2. 

BEN OUTPUT: DISTANCES DA A DB. AND AZIMUTHS ZA A ZB. 

GA^TANCPA-KSl/Z) :GB=TAN(P4-KS2/2) :DL=FNL(L1-L2) :DK»ABS(DL/BD) 
GA-GA-EP*(GA-0) :GB-GB-EP*(GB-0) 

DG>(G2-Gi)/BD:T=OL:X=LOG(GB)-LOG(GA) :GOSUB 160:ZA»A 
CZ«C08(ZA):DA-60«DE*C08(Gl):ir AB 8 (a)>. 0001 THEN 0A-60«DG/CZ 
E1-BC«SIN(G1) :E2-EC«SIN(G2) :ED-BC/2 

>DL:X-LOG(GB/( (l'»E2) /(i>E2)) *ED) -LOGCGA/ ((ItBi) / (1-El) ) *ED) 

GOSUB 160:ZB-A 

CZ-C08(ZB):DB-60«DK*C08(G1):IF AB8(CZ)>.0001THEN DB-60*DG/GZ 

BETUBN 

BEN 

BEN DIBECT SOLUTION. BHUMB LINE FOB SPHEBIGAL A SPHEBOIDAL EABTH. 

BEN INPUT: UTITUDE GI. LONGITUDE LI. FOBWABD AZIMUTH Ai AND 
BEN DISTANCE DO TO A POINT P2. NOTE: DD IS IN BADIANS. 

BDf OUTPUT: UTITUDE G2 AND LONGITUDE U A LS 
C6-C0S(A1) 

IF Ce^THEN LB^L1«DD/C08(G1):L8«LB:G2=G1:BETUBN 
G2-Gi«D0*C6 

GA«TAN(P4-^l/2) :GB«TAN(P4<KS2/2) :T6-TANUl) 

LB^L1-T6* (LOG(GB) -LOG(GA) ) 

Ei-BC«8IN(G1) :E>BC«8IN(G2) :ED-EC/2 

LS-Ll-ie* (LOG(GB/ ( ( 1+E2) /( 1-E2) ) -ED) -L0G(GA/( ( 1+El )/ (1-El) ) *ED) ) 

BETUBN 

REN 

REN DECIMAL TO DOD NN.F 

“:IF KOTHEN V$----:X— X 

X-X+1/1200:Y-INI(X) :V|-V|+RIGHT|(" -♦8TR$(Y) .3)+“®“ 
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1230 

1240 

1250 

1260 

1270 

1280 

1200 

1300 

1310 

1320 

2000 

2010 

2020 

2030 

2040 

2060 

2060 

2070 

2080 

2000 

3000 

3010 

3020 

3030 

3040 

3060 

3060 

3070 

3080 

3000 

3100 

3110 

3120 

3130 

3136 

3140 

3160 

3160 

3170 

3180 

3100 

3200 

3210 

3220 

3230 

3240 

3600 

3610 

3620 

3630 

3640 

3660 

3660 



»i600*(X-Y) tY^INTa) :H-8IE8(1000+Y) 

V8-V$+Mmai,3,2)+- . •+RICHT|(X|, 1)+*^' -.aETURN 
8Bi 

BBi ODD.NHSS TO DECIMAL 

nW:POI >1T0 LEM(V$):C|»MID$(V$,Z,1):IF 0-.-THE6 H-Z 
HSn :IF U-OTHEM Z-VAL(V<) :HET08N 
X-VALaKPT|(V|,IX)):81l-l:IF X<OTHEM 8N— 8K;X— X 
Vt»Vt+*0000* :T-VAL0IID$(V8. n+1 , 2)) :>VAL(MID8(V|, IX+3, 2)) 
X-8K*((Z/60+Y)/60^X) :BETURN 

CL8 :PEHfT 8PC(16); ‘NAVIGATION ALGORITHM 0E310 

PRINT :PRINT :PRINT 

PRINT SPG(16):‘l) DIRECT SOLUTION 

PRINT SPG(16);‘2) INVERSE SOLUTION 

PRINT SPG(16);‘3) RHOMB LINE APPROXIMATIONS" 

PRINT SPCC16) ; ‘ TO GREAT CIRCLE ROUTES 
PRINT :PRINT SPC(16);‘4) QUIT 

GOSUB 0010:G»VAL(C$):ON CG08UB 3000.3600,4000.6500 

GOTO 2000 

REN 

CL8 :PRINT SPC ( 16 ); ‘DIRECT SOLUTION": PRINT :PRINT 
PRINT SPC(IO);: PRINT ‘1ST UTITUDE DD.MN88 (-S) 

INPUT V|:GOSUB 1270:G1-X«RD 

PRINT SPC(IO);: PRINT ‘1ST LONGITUDE DDD.NM88 (-E) 

INPUT V8:G0SUB 1270:L1-X«RD 

PRINT SPC(IO);: PRINT ‘INITIAL COURSE DDD.MM88 

INPUT VI: GOSUB 1270:A1-X*RD 

PRINT SPC(IO); ‘DISTANCE (N.MI.) ";:INPOT D 

01-D*RD/60:DI><D/AE 

GOSUB 240: PRINT : PRINT SPC(8) ; ‘SPHEROID EARTH DIRECT SOLUTION 
NO-lOO: GOSUB 3210 

D1>-D1:G08UB 640:PRINT :PRINT SPC(8); ‘SPHERICAL EARTH DIRECT SOLUTION 
GOSUB 3210 

GOSUB 1040:PRINT :PRINT SPC (8) ; ‘RHUMB LINE SOLUTIONS"; 

PRINT SPC(4);‘UT‘;SPC(14);‘L0NG 
PRINT SPC(12); ‘SPHERICAL 

X-PNG(G2):X-X/RD: GOSUB 1210 :VS|-V|: PRINT Y|;SPC(8); 

X-FNL(LR):X-X/RO: GOSUB 1210:PRINT V| 

PRINT SPC(12) : ‘SPHEROID ‘;VS|;8PC(8); 

X-FNL(LS):X-X/RO: GOSUB 1210 : PRINT V| 

GOSUB 6000:GOTO 2000 
REN 

PRINT SPC(12);‘2ND UTITUDE :X»G2/RD:G0SUB 1210:PRINT VI 
PRINT SPC(12};‘2ND LONGITUDE :X-L2/RD:G0SUB 1210:PRINT V| 

PRINT SPC(12) ; ‘BACX AZIMUTH :X«A2/RD: GOSUB 1210:PRINT V|:RETURN 

REM 

CLS :PRINT SPC ( 16) ;‘ INVERSE SOLUTION": PRINT :PRINT 
PRINT SPC(IO);: PRINT ‘1ST UTITUDE DD.MNSS (-S) 

INPUT VI: GOSUB 1270:G1-X 

PRINT SPC(IO);: PRINT ‘1ST LONGITUDE DDD.NMSS (-E) 

INPUT VI: GOSUB 1270:L1-X 

PRINT SPC(IO);: PRINT ‘2ND UTITUDE DD.NMSS (<S) 

INPUT VI: GOSUB 1270 :G>X 
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3S70 

3680 

3690 

3000 

3010 

3020 

3030 

3040 

3060 

3000 

3070 

3080 

3090 

3700 

3710 

3720 

3730 

3740 

3760 

4000 

4010 

4020 

4030 

4040 

4060 

4000 

4070 

4080 

4090 

4100 

4110 

4120 

4130 

4140 

4160 

4100 

4170 

4180 

4190 

4200 

4210 

4220 

4230 

4240 

4260 

4200 

4270 

4280 

4290 

4300 

4310 

4320 

4330 



PRINT SPG(10);:P&INT •TSO LONGITUDE DDD.MM88 (-E) 

INPUT V$:G08UB 1270:L>X 
Gl-Gl*RD:G2-G2*RD:Ll-Ll*aD:L2-L2*BD 

GOSUB 440: PRINT : PRINT SPC(8) ; "SPHEROID EARTH INVEB8E SOLUTION 
U0-I00:n-1:G08UB 3700 

DD-D1:G0SUB 740:PRINT :PRINT SPC (8) SPHERICAL INVERSE SOLUTION 
IX=2: GOSUB 3700:G08UB 830 

PRINT : PRINT SPC(8) ; "BHUMB LINE SOLUTIONS: SPHERE SPHEROID 
PRINT SPGC12) ; "DISTANCE ";FNR(DA);SPC(8);FNR(DB);" N.NI. 

PRINT SPCC12) ; "COURSE "; :Z-ZA/RD: GOSUB 1210:PRINT V$;SPG(0); 
X«ZB/RD:GOSUB 1210:PRINT V< 

GOSUB 6000: GOTO 2000 
REM 

PRINT 8PCC12) ; "DISTANCE "; 

IF IX^ITHEM PRINT FNR(DD*AE);" N.NI. 

IF H-2THB1 PRINT FNR(60«DD/RD) ; " N.NI. 

PRINT SPCC12); "FORWARD COURSE "; :X-A1/BD: GOSUB 1210:PRINT V$ 

PRINT SPC(12);"BACX COURSE "; :X-A2/RD: GOSUB 1210:PRINT Y|:RETURN 
REN 

CLS : PRINT SPC(16); "RHUMB LINE APPROXIMATIONS": PRINT : PRINT 
PRINT SPC( 10);: PRINT "INITIAL UTITUDE DD.NM88 (-8) "; 

INPUT V$:G08UB 1270:G1-X«RD:G6-G1 

PRINT SPC(IO) ;: PRINT "INITIAL LONGITUDE DDD.NNS8 (-E) "; 

INPUT V|:GOSUB 1270:L1-X*RD:L6«L1 

PRINT SPC(IO);: PRINT "FINAL UTITUDE DD.NN88 (-S) "; 

INPUT V$:G08UB 1270:G2-X*RD:G7-G2 

PRINT SPC(IO);: PRINT "FINAL LONGITUDE DDD.NN88 (-E) "; 

INPUT VI.GOSUB 1270:L>X*RD:L7-L2 

GOSUB 740 :N0»10: PRINT .PRINT SPC(IO); "GREAT CIRCLE SOLUTION 
PRINT SPC(16); "INITIAL COURSE "; :AZ-A1:X-A1/RD: GOSUB 1210:PRINT V| 
PRINT SPCC16) : "DISTANCE ":FNR(eO*DD/RD);" N.NI. 

GOSUB 830:PRINT :PRINT SPC( 12) ; "RHUMB LINE (MERCATOR) SOLUTION 
PRINT SPCC16); "COURSE "; :X-ZB/RD: GOSUB 1210:PRINT VI 
PRINT SPCC16) ; "DISTANCE ";FNR(DB);" N.NI. 

GOSUB 4610 

PRINT : PRINT SPC(12); "VERTEX: UTITUDE "; :X-GV/RD: GOSUB 1210: PRINT V| 
PRINT SPCC20); "LONGITUDE "; :X-LV/RD: GOSUB 1210: PRINT V| 

GOSUB 6000 

CLS :PRINT SPC( 16) ; "RHUMB LINE APPROXIMATIONS" : PRINT 
PRINT : PRINT SPC(IO); "INPUT THE NAXINON NUMBER OF DEGREES 
PRINT SPC(IO);" OF LONGITUDE ON EACH RHUMB LINE LEG "; 

INPUT VI: GOSUB 1270:IN-X*RD 

G8-0:RS-0:IF IN-OTHEN IN-TP 

0L-FNL(L7-L6):IF OL-OTHEN F-1 

DV-FNL(LV-L6):IF DLoOTHEN F-DV/DL 

LX-0:IF F<-OOR F>-1THDI LX-1:G0T0 4370 

PRINT : PRINT SPC( 10) ; "LIMITING UHTUDE DD.MNS8 (-S) 

PRINT SPC(IO);" (0 TO OMIT) ";:INPUT V|: GOSUB 1270:LL-X*BD 
L8-0: L9-0 : G1-G6 :G3-G1 :L1-L6 : Ul-Ll : PH-LL : G2-G7 : L2-L7 
IF LL-OTHEM LX-1:G0TQ 4370 

qi-AB8(PV-G6):q>AB8(PV-G7):A-qi:IF q2<A THEN A-q2 
A-A«8GN(PV):B-PV-LL:IF A-0 THEN G-6:G0T0 4340 
G-B/A 
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4340 IF 6<-0 OR 6>1 THEN LW:PR1NT 8PC(10);*GANT USE. - IGNORED.*: GOTO 4370 
4360 G08DB 4810: PRINT 8FC(10); "LIMITING LONGITUDES: 

4360 Z"L8/SD:G08UB 1210:PRINT SPC(13) ;V|;SPG(10) :>L0/RD:G08UB 1210:PRINT V| 
4370 G08UB 6000:G08UB 5220 

4380 D8-FHL(L8-L1):D>FNL(L0-L1):IP AB8(D0)<AB8(D8)THEN T-L0:L>L8:L8-T 
4390 NX-1:LZ-L6:PX-G6:IF LX-ITHEN 4460 
4400 LT«L8:PT-LL:G08UB 4910 

4410 Ll«L8:61aLL;L2»L9:G2«LL:G0SUB 740:G0SUB 830 

4420 GOSUB 6110:GS-G8’»DD:RS-R8-»DB 

4430 LZ«L9:PZ-LL:L1-LZ:G1-FZ 

4440 L2>L7:G2«07 

4460 MX-0:LY-L7:PYNS7: GOSUB 4910 

4460 REM 

4470 PRINT : PRINT 8PC(10);*1) NEW APPROXIMATION OR 2) MEN PROBLEM? 

4480 GOSUB 6010:C»VAL(C|) :0N CG08UB 4190,2000 
4490 GOTO 4480 
4500 RSI 

4600 REM COMPUTE VERTEX 

4610 SA-aiN(AZ) :GA<08(AZ) :8P-8IN(G1) :CP-C0S(G1) 

4620 IF SA-OTHBl GV-P2:LV«L6: RETURN 

4630 LV*iU-ATN(l/8P/TAN(A2)):U-Ll:LV-FNL(LV):LM-LV:G08UB 4720 
4640 GV=>PH: RETURN 
4660 REN 

4700 RSI FIND LATITUDE GIVEN LONGITUDE 

4710 SA-SIN(AZ) :CA-COS(AZ) :SP-8IN(PA) :CP-C08(PA) 

4720 DL-LA-LM:Y-8P«C08(DL)*8A+SIN(DL)*CA:X-CP*8A: GOSUB 130:PH-A 
4730 PH-FNG(PR): RETURN 
4740 REM 

4800 REM FIND LONGITUDE GIVEN UHTUDE 

4810 n-0:IF AB8(PH)->AB8(GV)THEN KX-1:RETURN 

4820 SA-8IN(AZ):CA<08(AZ):SP-SIN(G3):CP-C08(G3) 

4830 E-CP*8A*TAN(PR):O8P*8A:S-CA:Ra-8qR(8*8tC*C):Y-8:X<:G0SUB 130:NU-A 
4840 COFNC(K/RO) :DD«LA1-NU:L8«FNL(00*CC) : L>FNL(D0-K:C) .RETURN 
4860 REM 

4900 BEN RHUMB APPROXIMATIONS SUBROUTINE. INPUT: LX. PX, LY A FT. 

4910 OL=FNL(LY-LZ) : IC=AB8(IN) :AG=AB8(DL)/IC 
4920 IF IC/RD<1THEN NP-INTaG>.5):M0«. 1 
4930 IF IC/RD>«1THEN NP-INT(AG) :N0»1 
4940 IF NP<1THEN IODL:L2»LY:GOTO 4970 
4950 IC-8GN(DL)*IC 

4960 L2-FNR(CLX+(DL-(HP-l)*IC)/2)/RD+.5*8GN(DL))*BD 

4970 L1»LX:G1»PX:LM»L2: GOSUB 4720:G2»PH 

4980 GOSUB 740:G0SUB 830:GOSUB 6110;GS»GS+DD:RS>RS'»DB 

4990 IF NP<1THEN 6060 

6000 N01:IF NP-ITHEN 6040 

6010 Ll-L2:Gl-G2:L2-FNL(Lli>IC):LM-L2: GOSUB 4720:G>PH 
6020 GOSUB 740:G0SUB 830:G0SUB 6110:G8=GS+D0:RS=RS-»DB 
6030 NC«NCi>l:IF NC<NPTHEN 6010 
6040 L1-L2:G1«G2:L2-LY:G2-PY 

6060 GOSUB 740: GOSUB 830; GOSUB 6110:G8-<»i-D0:R8-RSfDB:IF Mt-ITHEN RETURN 
6060 L1-LY:01-PY:M0-TP:A1-FNM(A2+PI): GOSUB 6110: RETURN 
6070 REM 

6100 BEN PRINT LINE OF OUTPUT 
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6110 

5120 

6130 

6140 

6160 

6160 

6170 

6180 

6100 

6200 

5210 

6220 

6230 

6240 

6260 

6260 

6270 

6600 

6610 

6000 

6010 

6020 

6030 



MO-1:X-01/BO:0080B 1210:PmT V$;SPC(2); 

X-L1/R0:G08UB 1210:PmT Y|; 

X-FNB(A1/BD):00SUB 6200:PB1NT V|;SPC(4); 

XaFNB(60«G8/B0) :G08UB 5200:PKINT V$;SPC(4): 

X-FNR(ZB/BO):G080B 6200:PSINT V$;SPC(3); 

X»FNR(R8) :G08UB 6200:PRINT V$ 

NL>NL-M:IF NL»20THEN G08UB 6000:GOSDB 6220 

RETUHM 

REN 

V$-* ■♦8TR$(X):V|-RICHT$(V$, 7) -.RETURN 
REM 

CIS :PR1NT SPCao); "RHUMB LINE APPROXINATION TO GREAT CIRCLE COURSE* 
PRINT 

PRINT * GREAT GREAT RHUMB RHUMB 

PRINT * CIRCLE CIRCLE LINE LINE 

PRINT * UT LONG COURSE DIST COURSE DIST 

NL-6:RETURN 

REM 

CLS .END 
REM 

PRINT :PRINT 8PC(10); "PRESS ANY KET TO CONTINUE 
FOR II-ITO 0:C$-INKEY$ iNEXT 
CS-INXEYR :IF C|-""THEN 6020 
RETURN 
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D* I^ogrflxn AmiotAtio]i^~NA!NrAXiCrOIt4 



Lme(s) 



Usage 



40 

50 



60 

70 

80 

90 

100 

110 

130 

140 

160 

240-350 

440-580 

640-680 

740-780 

830-900 



1040-1100 



1200-1240 

1270-1310 

2000-2060 

2070 

3000-3230 

3080 

3090 

3100 



PI 4 = f/4. PK = sr/2. PI = *•. TP = 2x. RD = t/ 180 is the degree-to-radian 
conversion. 

FL is the earth’s flattening factor. A£ is the earth’s equatorial radius. (Both 
of these constants are from the WGS-72 earth model.) ECC is the e^h’s 
eccentricity. 

FNM is the x mod MO function. 

FNL adjusts longitude to lie between -ir and k, 

FNR rounds x to the nearest MOf h. 

FNG adjusts latitudes to lie between -ir /2 and ir/2, 

FNACS Is the arccosine function. In the Sharp PC-1500A, function calls may 
be replaced inline by the ACS function. 

FNASN Is the arcsine function. In the Sharp PC-1500A, function calls may be 
replaced inline by the ASN function. 

FNATN 2 is the qatn function which returns a principle value between -ir and 

It, 

FNATNP is a qatn function which returns a principle value between 0 and 2x. 

GOTO 2000 to commence execution. 

Computation of the spheroid earth direct solution. Gl = Ll = Ai, A 1 =s 
aiQ and distance DD are input. G 2 = (^, L 2 = A^ and A 2 = a<ii are output. 

Computation of the spheroid earth inverse solution. Gl =s (^i , Ll = Ai , G 2 = 
^ and L 2 = Aq are input. A 1 = ai 3 , A 2 = a^i and distance DD are output. 

Computation of the spherical earth direct solution. Gl = LI = Ai, Al = 
ai 2 and distance DD are input. G 2 = L 2 = A^ and A 2 = a^i are output. 

Computation of the spherical earth inverse solution. Gl = Ll = Ai, G 2 = 
<fh 2 and L 2 = A^ are input. Al = ais, A 2 = a^i and distance DD are output. 

Computation of the rhumb line inverse solution. Gl = Ll = Ai, G2 = 
and L 2 = A^ are mput. Spherical earth distance DA and asimuth ZA as well 
as spheroid earth distance DB and asimuth ZB are output. 

Computation of the rhumb line direct solution. Gl = Ll = Ai, Al = 
ati 3 and distance DD are mput. G 2 = ^ 3 , spherical earth longitude LR and 
spheroid earth longitude LS are output. 

Convert decimal degrees to degrees, minutes and tenths of minute format. 

Convert packed degrees, minutes and seconds format (DDD.MMSS) to decimal 
degrees. 

Master menu option display. 

lYansfer to selected menu item. 

Direct solution option. 

D is the distance (input) in nautical miles, D1 is the spherical earth distance 
in radians and DD is the distance in earth equatorial radius units. 

Compute spheroid earth direct solution. Print heading. 

Print spheroid earth direct solution. 
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3110 

3120 

3130 

3140-3160 

3170-3180 

3190 

3210-3230 

3500^740 

3600 

3610 

3620 

3630 

3640 

3650 

3660-3670 

3680 

3700-3740 

4000-5260 

4000-4080 

4090-4110 

4120-4140 

4150-4170 

4180 

4190-4220 

4230 



4240-4260 

4270-4280 

4290-4340 



4350-4360 

4370 

4380 

4390 



Compute spherical earth direct solution. Print heading. 

Print spherical earth direct solution. 

Compute rhumb line solution. Print heading. 

Print spherical earth rhumb line latitude and longitude. 

Print spheroid earth rhumb line latitude and longitude. 

Continue prompt. Go to master menu. 

Output subroutine for option 1. 

Inverse solution option. 

Compute spheroid earth inverse solution. Print heading. 

Print spheroid earth inverse solution. 

Compute spherical earth inverse solution. Print heading. 

Print spherical earth inverse solution. 

Heading for rhumb line solutions. 

Print spherical and spheroid earth rimmb line distances. 

Print spherical and spheroid earth rhumb line courses. 

Continue prompt. Go to master menu. 

Output subroutine for option 2. 

Rhumb line approximation to great circle option. 

hiput prompts. 

Compute great circle solution and print initial course and distance. 

Compute rhumb line solution and print course and distance. 

Compute and print the latitude and longitude of the vertex of the great circle 
route. 

Display continue prompt. 

Prompt for the longitude increment of the rhumb line approximation. 

GS and RS are the cumulative great circle and rhumb line distances traveled on 
each leg, respectively. Initialize them to zero. IN is the longitude increment 
in radians. H input as zero (no increments requested), it is set to 2x. 

Determine if the vertex lies between the origin and destination. If not, go to 
4370. Otherwise proceed. 

H the vertex is on the great circle route, prompt for a limiting latitude LL. 

Set limiting longitudes L8 and L9 equal to zero. Determine if the limiting 
latitude cuts the great circle course. If it does not, inform the user that the 
limiting latitude is ignored. L% is zero if the limiting latitude is to be used, 
otherwise it is one. 

If the limiting latitude cuts the great circle course, compute and display the 
longitudes L8 and L9 at which the limiting latitude cuts the great circle. 

Display continue prompt, then print heading. 

Determine which limiting longitude is closest to the initial position. L8 becomes 
the closest, L9 the farthest. 

LX and PX are the initial great circle longitude and latitude. If the limiting 
latitude b not to be used, go to 4450. 
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4400 

4410 

4420 

4430-4440 

4450 

4470-4490 

4610-4640 

4710-4730 

4810-4840 

4910-5060 

4910-4950 



4960-4980 

4990 

5000 

5010-5030 

5040-5050 

5060 

5110-5200 

5220-5260 

5500 

6000-6030 



LY and PY are the longitude and latitude of the final point on the first segment 
of the great circle coarse. GOSUB 4910 to compute the rhumb line apprcudma- 
tion up to the first limiting longitude. 

Compute the great circle and rhumb line distances from the limiting latitude 
at the first limiting longitude to the second limiting longitude. 

Print one line of output (for the limiting longitudes). Update the distance 
counters. 

Reset LX, PX, LY and PY for the section of the great circle and rhumb line 
following the limiting latitude section of the course. 

GOSUB 4910 to compute the final section of a course following a limiting latitude 
leg or compute the entire course for cases in which the limiting latitude is not 
used. 

Prompt for either a rework of the rhumb line approximation with new param- 
eters or for a new problem. 

Compute the latitude GY an longitude LY of a great circle vertex. 

For a given longitude, compute the corresponding latitude of on a great circle. 

For a given latitude, compute the corresponding longitudes L8 and L9 on a 
given great circle. K% is set to one if there is no solution, otherwise is 
zero. 

This is the subroutine which selects the longitude endpoints of the rhumb line 
approximation to a great circle between a starting position with latitude PX 
and longitude LX and a final position with latitude PY and longitude LY. It 
then computes and prints the heading and cumulative distance along each leg. 

Using the increment size, IN and IC, compute the number of interior points, 
NP, on the rhumb line approximation exclusive of the initial and final points 
(the number of interior legs is NP — 1. Redefine IC if it is longer than the 
difference of the initial and final longitudes. 

Find the heading and distance on the initial leg. Print out first line of output. 

If there are no internal legs, go to 5060. 

Initialize the counter IC. K there is only one leg remaining, go to 5040. 

Loop to compute print each internal leg. 

Compute the headings and distances for the final leg. If M% is one, return to 
compute the limiting latitude leg. 

Print the final summary line of output. 

Subroutine to print one line of output. NL counts the number of lines of output. 
If 20 lines are output, a new screen is started. 

Subroutine to print screen heading. 

Return to BASIC. 

“PRESS ANY KEY TO CONTINUE^’ subroutine. 
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IV. NAVEPHMt Almaiuic and Ephemeria Program. 

A. Introduction. The basis for NAVEPHM are the equations of VanFlandem and Pulkki- 
nen (VFP) [Ref. 10] . These equations can be used to compute the mean heliocentric 
positions of the sun and planets and the mean geocentric position of the moon for the 
mean ecliptic and equinox of date. The authors claim their formulas to be of low-precision 
(1^) and valid for any epoch within 300 years of the present. When corrected for nutation 
and aberration the accuracy of their formulas, at least for the sun, moon and navigational 
planets, appears to be much better, i.e, 0.2^, at least. 

The formulas of VFP are used to compute heliocentric spherical ecliptic coordinates 
for any specified q}hemeris time, ET. These coordinates are longitude A, latitude ^ and 
radius vector r. These coordinates must be converted to rectangular coordinates x, y, and 
z using the standard transformation 

X SB r cos COS A, 

y = fcosj^sinA, and (30) 

To obt^ the geocentric x, y, and z coordinates of the planets, subtract the x, y, and z 
coordinates of the sun from the x, y, and z coordinates of the planets, respectively. 

The Julian day number, required in many calculations, is obtained using the equation 
on page B2 of Ref. 2. The ephemeris time is obt^ed from the universal time, UT, from 
the equation £T s UT + AT. The factor AT is obtained by astronomical observation 
only. The formula used here, AT = 81.94T - 15, was obtained by least squares from the 
1980.5 through 1985.5 values given on pages B5 and K9 of The Astronomical Almanac 
1985 [Ref. 11], where T is the number of Julian centuries elapsed frrom 1900 January 0, 
12^ ET. Although a plot of AT as a function of time is linear for 1980.5 through 1985.5, 
this should be checked with each new edition of The Astronomical Almanac. An accurate 
value of AT affects only the computation of the moon’s position. Errors of as much as 
9.6* in AT will affect the computation the moon’s position by at most 0.1’, consequently 
it will not be necessary to change FMPL distribution tapes yearly. 

The heliocentric in- plane velocity components, x<j and y^, of the planets, required for 
the aberration correction, can be computed from the formulas [Ref. 12, pg. 85]: 

Xu, = —oeEsmE, and 
y„ = aE\/l - e^cosE, 
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where E is the eccentric anomaly, ^ is the time derivative oiEya ia the semimajor axis and 
e is the eccentricity. Unfortnnately, to find E one mnst solve Kepler’s equation iteratively, 
which is a slow process. With an error of no more than 2% for the navigational planets, 
one can use the mean anomaly M in place of ^ in Equs. (31). The in-plane, velocity 
components are thus approximated by 

= -oeE^My and . . 

. , (32) 

= aEy/l — e^cosM. 



IVom Equ. (3.76) of Ref. 12, it can be shown that 







where /i is the reduced mass and k is the gravitational constant. The heliocentric ecliptic 
velocity components, x, y and i, can then be obtained by the transformation 
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where Q is the longitude of the ascending node, i is the orbital inclination and u/ is the 
argument of perihelion. The equations for computing, a, e, Af, a;, i and 0 were obtained 
from Escobal [Ref. 13, pp. 8-9|. The geocentric i, y and z coordinates of the planets are 
obtained by subtracting the z, y and i coordinates of the sun from the heliocentric z, y 
and z coordinates of the planets. Since the aberration of the moon is negligible, its velocity 
components are not computed. 

Once the mean geocentric rectangular positions (zm, thn ^d Zm) ^d velocities (zmi 
Hm and im) have been obtained, the longitude and latitude for the mean ecliptic and 
equinox of date can be obtuned by inverting Equs. (30), so that 



®m)i and 
= sin-* = tan-‘ 



2m 






(33) 



Next the nutation of longitude and obliquity Ae are computed [Ref. 4, §2C]. The 
geocentric rectangular positions (zt, yt and zt) for the true ecliptic and equinox of date 
can be obtained from the transformation 
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The same transfonnation la used for the velocity components. The longitude and latitude 
for the true ecliptic and equinox of date can be obtained by substituting the true positions 
into Equs. (33). 

The aberration or light-time correction converts true positions into apparent positions. 
The x-coordinate transformation is 







where r is the geocentric distance and c is the velocity of light. Similar transformations 
are used for the y- and z-coordinates. The apparent longitude and latitude for the true 
ecliptic and equinox of date can be obtained by substituting the apparent positions into 
Equs. (33). 

The apparent geocentric rectangular positions (xaqi Vaq ^^d Zaq) for the true equator 
and equinox of date are obtained from the transfonnation 
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The apparent right ascension a and declination 6 for the true equator and equinox of date 
can be computed from 



a = qatn(y<i9, £04)) ^d 

« = ain-‘ = tan-* , ^ («) 

*■ + 

The Greenwich mean sidereal time, Greenwich apparent sidereal time, Greenwich 
hour angle and local hour angle are computed from formulas given on pages B3 and B4 
of Ref. 2. The altitude and azimuth are computed from Equs. (10) and (11), respectively. 
The equation for refraction is Equ. (3), page B15 of Ref. 2. The equations for planetary 
magnitude are given on pg. 315 of Ref. 4 with the correction for Saturn’s rings given on 
pages 362 to 365. Formulas for the table on page 365 were obtained using least squares. 
Equations for the semidiameter of the sun and moon are found on page B16 of Ref. 2. The 
lunar parallax in altitude formula is on page B16 of Ref. 2 and the lunar phase formula is 
on page 311 of Ref. 4. The lunar age approximation developed by the author is within ±1 
day. 
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B. Sample Problem. Calculate the position of the sun, moon and navigational plan- 
ets for 24 May 1984 at 17^58^45^ zulu time at latitude 38*^35^24" north and longitude 
76^18W' west. Assume that the temperature and pressure are the default values of 10^ C. 
and 1010 mb., respectively. The output shown in this sample problem was computed using 
Microsoft BASIOA/D (double precision) on an IBM PO. Several intermediate results, such 
as the true position, and apparent positions for the true equator and for the true equinox 
of date are printed only as a debugging md for those wishing to reproduce the results on 
a computer other than those for which NAVEPHM is available. 

Input Screen: 

YEAR (4 DIGITS) 7 1984 

MONTH NUMBER 7 5 

DAY OF THE MONTH 7 24 
ZULU TIME (HH.MMSS) 7 17.5846 
LAT (DD.MMSS) (4N/-S)7 38.3624 
LON (DDD.MMSS) (•>■¥/-£) 7 76.18 
TEMP (DEG CELSIUS) 7 10 
PRESSURE (MILLIBARS) 7 1010 

PRESS C TO CONTINUE 

The input screen is used for input only — there are no options from which to choose. The 
temperature and pressure are used only for the atmospheric refraction computation. 
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1st Ontpni Screen: 

SUN 

TRUE POS, TRUE ECL k EQNX OF DATE (UT) : 

63.73236366163271 0 

APP POS. TRUE ECL k EQNX OF DATE (UT) : 

63.72674638386791 0 

APP POS. TRUE EQU k EQNX OF DATE (UT) : 

4 . 11446999613666 20 . 89927908414662 

GREENWICH HOUR ANGLE k DECL (UT) : 

90^29.0' 20®64.0^ 

ALTITUDE - 68.60634784607299 68^30.3' 

AZIMUTH » 218.6611413863376 218<’39.r 

SD - 16.8^ 

REFRACTION - .4' 

PRESS C TO CONTINUE 

The trne and apparent poations for the true ecliptic and equinox of date printed are the 
ecliptic longitude and latitude, respectively, in decimal degrees. The apparent position 
for the true equator and equinox of date is the right ascension in decimal hours and the 
decimation in decimal minutes. The Greenwich hour angle and declination are given in 
degrees, minutes, and tenths of minute notation. The altitude and asimuth are given in 
both decimal degrees and degrees, minutes, and tenths of minute notation. SD is the sun’s 
semi-diameter. The refraction is for the computed altitude of the sun. 
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2nd Ontpnt Screen: 
KOON 



TRUE PCS. inUE ECL k EQNX OF DATE (UT) : 

366 . 6933313649804 -4 . 991896604201667 

APP POS. TRUE ECL k EQNX OF DATE (UT): 

366 . 6933313649804 -4 . 991896604201667 

APP POS. TIUJE EQU k EQNX OF DATE (UT) : 

23 . 86924633629027 -6 . 291911880374428 

GREENWICH HOUR ANGLE k DSa (UT) : 

164*09.7^ - 6®17.6^ 

ALTITUDE - 6.46172646204884 6®27.1' 

AZIMUTH > 267.4666216802323 267'’ 28.0' 

PHASE .31 AGE - 24 DAYS 
SD - 14.8' SD AUG - 14.8' 

P IN A - 64' 

REFRACTION - 10.2' 

PRESS C TO CONTINUE 

The output for the moon is similar to that for the sun. In addition, the moon's phase, 
approximate age, augmented semi-diameter (topocentric semi- diameter) and parallax in 
altitude P IN A are output. The true position and the apparent position for the true ecliptic 
and equinox of date for the moon are identical because no aberration correction is made 
for the moon (see Section A). 
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3rd Output Screen: 
VENUS 



TRUE PQS, TRUE Ea k EQNX OF DATE (UT) : 

57.72453662418941 -.6628498164408445 

APP POS. TRUE ECL k EQNX OF DATE (UT) : 

57 . 71230645494549 - . 6532096744515719 

APP POS, TRUE EQU k EQNX OF DATE (UT) : 

3 . 706666485866183 19 . 01586622168487 

GREENWICH HOUR ANGLE k DEa (UT) : 

oe^oe.o' 19*01.0' 

ALTITUDE - 63.67703088674409 63*40.6' 

AZIMUTH - 227.7063426077327 227*42.4' 

MAGNITUDE - -3.4 
REFRACTION - .6' 

PRESS C TO CONTINUE 

The ontpnt for Venus is sinular to that for the sun except that the apparent magnitude is 
output instead of the semi-diameter. 
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4tli Ontpnt Screen: 
MARS 



TRUE POS, TRUE ECL k EQNX OF DATE (UT) : 

226 . 1549997906666 - . 7543627638263081 

APP POS. TRUE ECL k EQNX OF DATE (UT) : 

226.1661434601826 - .7542106991416879 

APP POS. TRUE EQU k EQNX OF DATE (UT) : 

14 . 89744688277316 - 17 . 39602620686829 

GREENWICH HOUR ANGLE k DECL (UT) : 

288*44.3' - 17*23.8' 

ALTITUDE - -64,68434644133616 - 64*41.1' 

AZIMUTH - 62.30703743162624 62*18.4' 

MAGNITUDE - -1.6 
REFRACTION - -.7' 

PRESS C TO CONTINUE 

The output for Mars is similar to that for Venus. Note that the altitude of Mars is negative, 
that is, Mars is below the horison. The leads to the anomolous computation of a negative 
value for the refraction. When a planet is below the horizon, the value computed for 
refraction, is meaningless. 
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Sth Ontpnt Screent 

JUPITZR 

mjz POS, TIUJE ECL k EQNX OF DATE (UT) : 

282.0081844049666 .1412184643412081 

APP POS, TRUE ECL k EQNX OF DATE (UT) : 

282 . 0100408679675 . 1412662061186686 

APP POS. TRUE EQU k EQNX OF DATE (UT): 

18 . 86941225466823 >22 . 76882790181082 

GREENWICH HOUR ANGLE k DEa (UT) : 

229*09.6^ - 22®46.6' 

ALTITUDE - -61.97026716907849 - 61® 68. 2' 

AZIMUTH - 296.4711031461716 296®28.3' 

MAGNITUDE - -2.1 
REFRACTION - -.6' 

PRESS C TO CONTINUE 

The output for Jupiter is similar to that for Venus. See the comments for Mars regarding 
a negative refraction value. 
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0th Ontpvt ScTe€n: 

SATURN 

TRUE PCS. TRUE ECL A EQNX OF DATE (UT) : 

221 . 6236884013894 2 . 680801361264901 

APP POS, TRUE ECL A EQNX OF DATE (UT) : 

221 . 6271666166949 2 . 68091611296868 

APP POS, TRUE EQU A EQNX OF DATE (UT) : 

14 . 66067626686222 - 12 . 83630963737889 

GREENWICH HOUR ANGLE A DEa (UT) : 

292®17.4^ - 12®60.2^ 

ALTITUDE - -49.04279937791162 - 49*»02.6' 

AZIMUTH » 60.93733261413396 60" 56. 2' 

MAGNITUDE - .4 
REFRACTION - -.8' 

PRESS C TO CONTINUE 

The output for Saturn is similar to that for Venus. See the comments for Mars regarding 
a negative refraction value. 
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C. Program Liating — NAVEPHM 



10 BQI YANFLANDERN A PULZKH1EN PLANET EPHD1ER18 . 07-20-86. 

REV 08-27-86 0 0800 

12 RDI "AOTRO EPHEN Y* 

13 REH •NAYEPHM.BAS* 

20 DEFCBU-N.O-Y:Dllf A(26).B$(6) 

30 P2-2*ATM(1) :PI-P2+P2:TP-PI+PI;RD-PI/180:8D-8D/3600 
40 EK«.01720200806:MB|«>1010*:DC$-‘‘10*:0N3-CHR|(26) 

60 DEPFNA(X)-(X-m(l))*TP:DEFPMB(I)-X-TP*IllT(X/TP) : 

DEFFNR(X)-INT(X*MO-» . 5)/N0 
60 DEFFNN(X)»X-M0*1NT(X/N0) 

70 B$(0)-«8UN« :B|(1)-**M00M" :B|(2)-"VENUS« :B$(3)-"MAR8« 

80 B|(4)-* JUPITER* :B6(5)«<8ATURN* 

90 GOTO 2090 
100 : 

no REH 2-ARG ATAN FCTN. RA01AN8. 

120 A>ATN(Y/ a- lE-09* (X-O) ) ) -PI* (X<0) +TP* (X>0) * (Y<0) : RETURN 
130 : 

140 REH ARC8IN A ARCC08 FCTN*8 
160 A8-ATN(X/(8QR(1-X*X)-1E-09*(AB8(X)-1))) :RETURN 
160 AOATN(8QR(l-X*X)/(X-lE-09*(X-0)))-PI*a<0) :RETURN 
170 : 

180 REM HELIO 8PHER-T0-RECT * NOTION COMPUTATION 
190 RpR/100000:C8-R*C08(B) 

200 XE-C8*C08(L) : YH-C8*8IN(L) : ZE-R*8IN(B) 

210 IF N-ITHEN RETURN 

220 MU-1: IF NOOTEEN NU-1*1/RN(N) 

230 PA-EX*8qR(AN(N)*NU)/R 

240 X— FA*8IN(M(N)) :Y-FA* 8 qR(l-BC(N)'' 2 )*C 08 (M(N)) :Z-0 
260 A— AP(H):GOSUB 280:A— IN(N) :GOSUB 290:A— AN(N) :GOSUB 280 
260 UH-X:YH-Y:yH-Z:RETURN 
270 : 

280 CA-C08(A):8A-8IN(A):>X*CA-»Y*SA:Y-Y*GA-X*SA:X-T:RETURN :REM Z-AXI8 ROT 
290 CA»C08(A):8A=8IN(A):T-Y*a-»Z*8A:Z=Z*GA-Y*8A:Y»T:R£TURN :REM X-AXIS ROT 
300 : 

310 REM RECT-TO-SPHER 
320 GOSUB 120 

330 L-A2/RD:R>X*X*Y*Y;B-ATN(Z/8qR(R2))/RD:R-8qR(R2<»Z*Z) :RETURN 
340 : 

360 REM DECIMAL TO ODD NM.F 
360 VI-* ■:!? X<OTHEN V$-"-":X— X 

370 X=X+1/1200:Y=INT(X):V6»V$+RIGHT$(* ■+8TR$(Y) ,3)+*®- 
380 X-600* (X-Y) : Y-INT(X) :X|-8TE8( 1000+Y) 

390 Y$»V|-»NID|(X8.3.2)-^*.*-^RIGHT|(X|.l):RETURN 
400 : 

410 REM DDD.NM8S TO DECIMAL 

420 IX-0:FOR Z-ITO LEN(V|):C$-MID|(V|,Z,1):IF C$-*.“THEH IX-Z 

430 NEXT :IF H-OTHEM X-VAL(V$) iRETURN 

440 X-VAL(LEFT$(V|.IX));SN-1:IF X<OTHEN 8N— 8N;X— X 

460 V8»V$+*0000* : Y»VAL(MID| (V| , IX+ 1 . 2) ) : Z»VAL(MID| (V| , IX+3 . 2) ) 

460 X-8N* ( (Z/e0+Y)/60+X) : RETURN 
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RDf FUNDAMENTAL AHGDNEKTS 

A(1)-.606434-».03660110120*TS:A(2)-.374807-^. 03620164700*78 
A(3)». 260001+ . 0367481062*78 :A(6)=0:A(7)=. 770072+ . 00273700031*78 
A(8)- . 003126+ . 0027377786*78 

A(12)«. 606408+4 .460468670000003E-03*78 :A(13)-. 140023+ . 00446036173*78 
A( 14)- . 202408+ . 00446040017*78 : A(16)- . 087363+ . 00146676328*78 
A( 16)- . 063866+ . 00146661327*78 : AC17) 840604+ . 00146660466*78 
A(18)- . 080608+ . 00023080803*78 :A(10)- . 066631+ . 00023080803*78 
A(20)-. 814704+ . 00023080803*78 :A(21)- . 133206+ . 00000204371*78 
A(22)- . 882087+ . 00000204371*78 : A(23)- . 821218+ . 00000204371*78 
A(26) -. 400680+3 . 260438000000001E-06*78 
A(4)-A(l)-A(7) :A(6)-A(1)-A(3) 

FOR Z-170 26:A(Z)-FNAa(Z)):NEZT Z:RE7URN 

raiNT B$(0):RO( SUN(O) 

X=A(8):Y=A(13):Z=A(10) 

L-(6010-17*7J)*8INa)+72*8IN(2*X)-7*C08(X-Z)+6*8IN(A(l)-A(7)) 
L-L+6*SIN(4*X-8*A( 16) +3*Z) -6*C08(2* (X-Y) ) -4*SIN(X-Y) 
L-L+4*C08(4*X-8*A( 16) +3*Z) +3*8IN(2* (X-Y) ) -3*SIN(Z) -3*8IN(2* (X-Z) ) 
L=L*8D+A(7) :B=0 

R^100014-1676*G08(X) -14*CQ8(2*X) 

AM(0)-1 . 00000023: EC(0)-. 016700114 :M(0)^(8) 

AP(0)-( (4 . 62778E-04*7J+1 . 710176) *7J+101 . 2208333) *RD: 

IN(0)-0:AN(0)-180*RD 

RN(0) -1/320300 

G08UB 100:RR-R:X8-XE:Y8-YH:Z8-ZH:U8-UH: :V8-VH:W8-«H: RETURN 
PRINT B|(1):RD( MQQN(l) 

W»A(2) :X=A(3) :Y=A(4) :Z»A(8) 

L-22640*8IN(W) -4686*SIN(W-2*Y) +2370*8IN(2*Y) +760*SIN(2*W) 

L-L-668*8IN(Z) -412*8IN(2*X) -212*8IN(2* (W-Y) ) -206*8IN0l»- 2*Y+Z) 

L-L+102*8IN(W+2*Y)+16B*8IN(2*Y-Z)+148*8IN(W-Z)-126*8IN(Y) 

L-L-110*8IN(W+Z)-66*8IN(2*(X-Y))-4B*8IN(W+2*X)+40*8IN(W-2*X) 

L-L-38*8IN(W-4*Y)+36*8IN(3*W)-31*8IN(2*W-4*Y)+28*8IN(W-2*Y-Z) 

L=L-24*8IN(2*Y+Z)+10*8IN(W-Y)+18*8IN(Y+Z)+16*8IN(W+2*Y-Z) 

L-L+14*8IN(2* (W+Y) )+14*8IN(4*Y) -13*8IN(3*W-2*Y) 

AG-W+16*A(7)-18*A(12):L»L-11*8IN(AC)+0*C08(AG)+4*TJ*(C08(AG)+8IN(AG)) 

L-L+10*8IN(2*W-Z)+0*8IN(W-2*X-2*Y)-0*8IN(2*W-2*Y+Z) 

L=L-8*8IN0f+Y)+8*8IN(2*Y-2*Z)-8*8IN(2*W+Z)-7*8IN(2*Z) 

L-L-7*8IN(W-2*Y+2*Z)+7*8IN(A(6))-6*8IN(W-2*X+2*Y) 

L-L-6*8IN(2*X+2*Y) -4*8IN(W-4*Y+Z) -4*8IN(2*W+2*X) 

L-L+3*(8IN(W-3*Y)-8IN(W+2*Y+Z)-8IN(2*W-4*Y+Z)+8IN(W-2*Z) 

+8IN(W-2*Y-2*Z)) 

L-L+2*(8IN(W+4*Y)-8IN(2*W-2*Y-Z)-8IN(2*X-2*Y+Z)) 

L-L+2*(8IN(4*W)+8IN(4*Y-Z)+8IN(2*W-Y)) 

B-18461*8IN(X) +1010*8IN(W+X)+1000*8IN(W-X) -624*8IN(X-2*Y) 

B=B-100*8IN(W-X-2*Y)-167*8IN(W+X-2*Y)+117*8IN(X+2*Y) 

>B+62*8IN(2*W+X)+33*8IN(W-X+2*Y)+32*8IN(2*W-X)-30*8IN<X-2*Y+Z) 

B-B-16*8IN(2*W+X-2*Y)+16*8IN(W+X+2*Y)+12*8IN(X-2*Y-Z) 

B-B-0*8IN(W-X-2*Y+Z)-8*8IN(X+A(6))+8*8IN(X+2*Y-Z) 

B-B+7*(-BIN(W+X-2*Y+Z)+8IN(W+X-Z)-8IN(W+X-4*Y)) 

B-B+6*(-8IN(X+Z)-8IN(3*X)+8IN(W-X-Z)) 

B-B+6*(-8IN(X+Y)-8IN(W+X+Z)-aiN<W-X+Z)+8IN(X-Z)+8IN(X-Y)) 
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>B+4* (8IN(3*W+X) -8INa-4*Y) ) +3* ( -SIN(W-X-4* Y) +Sm (W-3*X) ) 
>B+2*<-8IM(2*W-X-4*Y)-8IlI(3*X-2*Y)+8IN(2*W-X+2*Y)) 

B-B+2* (SIN(¥-X+2*Y-Z) +SIN (2*W-X-2*Y)+8IN (3*W-X) ) 

Bf6036208-327746*C08 (V) -670O4*C0S(W-2*Y) -46357*008 (2*Y) 

RfR-8004*CQS (2*W) *3865*008 (2*W-2*Y) -3237*008 (2*Y-Z) 
>R-2688*008(W+2*Y)-2358*008(W-2*Y+Z)-2030*008(W-Z) 

>8+1710*008 (Y) +1671*008 (W+Z) + 1247*008 (W-2*X) +704*008(Z) 

>R+620*008 (2*Y+Z) - 624*008 (W-4*Y) +308*008 (W-2*Y-Z) -366*008 (3*W) 
>8-205*008 (2*W-4*Y) - 263*008 (Y+Z) +240*008 (3*W-2* Y) -221 *008 (W+2*Y- Z) 
>8+185*008(2*X-2*Y)-161*OOa(2*Y-2*Z)+147*008(W+2*X-2*Y)-142*008(4*Y) 
>8+ 130*OOS(2*W-2*Y+Z) - 118*008 (W-4* Y+Z) - 116*008 (2*W+2*Y) 
>8-110*008(2*W-Z) 

>L*SD+A(1) :>B*8D: >8/23464. 8 :G0SUB 100:8P-8:8ETU8N 
PRINT B$(2);8EM VENU8(2) 

X-A(8):Y-A(13);>A(14) 

L=(2814-20*TJ)*aiN(Y)-181*8IN(2*Z)+12*8IN(2*Y)-10*008(2*(X-Y)) 

>L+7*008(3*(X-Y)) 

>1221S*8IN(Z)+83*(SIN(Y+Z)+8IN(Y-Z)) 

8^72336-403*008 (Y) 

>L*8D+A(12) :B*B*8D 

AN(2)- . 7233316000000001 : E0(2)-6 . 773041000000001E-03: N(2)-A(13) 
AP(2)»((- .001386380*TJ+. 60818611)*TJ+64.3841861)*8D: IN(2)»3.3046*RD 
AN(2)-( ( . 00041*TJ+ .80085) *TJ+75 . 77064722000001)*8D:8N(2)-408623. 6 
GOSOB 100:8P«8:BETI]8N 
PRINT B$(3) :8EU MA88(3) 

X-A(8) :Y-A(16) :Z-A(17) :Q-A(10) 

L-(38461+37*TJ)*8IN(Y)+(2238+4*TJ)*8IN(2*Y)+181*8IN(3*Y)-52*aiN(2*Z) 
L=L-22*008(Y-2*q)-10*8IN(Y-q)+17*008(T-q)+l7*8IN(4*Y)-16*008(2*(Y-q)) 
L»L+13*008(X-2*Y)-10*8IN(Y-2*Z)-10*8IN(Y+2*Z)+7*008(X-Y)-7*008(2*X-3*Y) 
>L-5.*8IN(A( 13) -3*Y) -S*8IN (X-Y) -5*8IN (X-2*Y) -4*008 (2*X-4*Y) +4*008 (Q) 
L-L+3*008 (A( 13) -3*Y) +3*8IN (2* (Y-q) ) 
B-6603*8IN(Z)+622*8IN(Y-Z)+616*8IN(Y+Z)+64*8IN(2*Y+Z) 

>153031 -14170*008(Y) -660*008 (2*Y) -47*008 (3*Y) 

>L*8D+A(15):>B*8D 

AN(3)-1 . 62368830: E0(3)>. 00340488700000002 :M(3)-A(16) 

AP(3)-( (1 . 3126E-04*TJ+1 . 06076667) *TJ+285 . 4317610000001 ) *BD : 
IN(3)-1.8407*8D 

AN(3)-( (-1 . 380E-06*TJ+ . 7700016700000001) *TJ+48 . 78644167) *RD : 

8M(3) >3008710 

GOSUB 100:RP>R:B£TDRN 

PRINT B|(4):8E3f JUPITE8(4) 

X-A(18) :Y-A(10) :Z-A(22) 

L-10034*8IN(Y)+5023*TJ+2611+1003*C08(2*Y-5*Z)+601*8IN(2*Y) 

L>L-470*8IN(2*Y-6*Z)-185*8IN(2*Y-2*Z)+137*8IN(3*Y-6*Z)-131*8IN(Y-2*Z) 

L-L+70*C08(Y-Z)-7e*C08(2*Y-2*Z)-74*TJ*C08(Y)+68*TJ*8IN(Y) 

+66*008 (2*Y-3*Z) 

L-L+63*C08(3*Y-6*Z)+63*C08(Y-6*Z)+40*8IN(2*Y-3*Z)-43*TJ*8IN(2*Y-6*Z) 

L=L-37*C08(Y)+25*8IN(2*X)+25*8IN(3*Y)-23*8IN(Y-6*Z)-10*TJ*C08(2*Y-5*Z) 

L>L+17*C08(2*Y-4*Z)+17*C08(3*Y-3*Z)-14*8IN(Y-Z)-13*8IN(3*Y-4*Z) 

-0*008 (2*X) 

L>L+0*008(Z)-0*8IN(Z)-0*8IN(3*Y-2*Z)+0*8IN(4*Y-5*Z)+0 

*8IN(2*Y-6*Z+3*A(25)) 

L>L-8*008(4*Y-10*Z)+7*008(3*Y-4*Z)-7*008(Y-3*Z)-7*SIN(4*Y-10*Z) 
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1460 L«L-7*8INCr-3*Z)+6*C08(4*Y-6*Z)-6*8IN(3*Y-3*Z)+5*C08(2*Z) 
-4*8IN(4*Y-4*Z) 

1470 L«L-4*CQ8(3*Z)+4*C08(2*Y-Z)-4*C08(3*Y-2*Z)-4*TJ*C08(2*Y)+3*TJ*8IN(2*Y) 
1480 L-L+3*C08(8*Z)+3*C08(5*Y-10*Z)+3*8IN(2*Z)-2*8IN(2*X-Y)+2*8I1I(2*X+Y) 
1400 >L-2*TJ*8IN (3*Y-6*Z) -2*TJ*8IN(Y-6*Z) 

1600 B— 4602*C08(Y)+2EO*8IN(Y)+227-227*C08(2*Y)+30*TJ*8INa)+21*TJ*C08(Y) 
1610 B-B+16*8IN(3*Y-6*Z)-13*8IN(Y-6*Z)-12*C08(3*Y)+12*8IN(2*Y) 

1620 B»B+7*C08(3*Y-6*Z)-6*C08(Y-6*Z) 

1630 Bp620883-26122*CQ8(Y)-604«C08(2«Y)-^260«CQ8(2«(Y-Z))-170«C08(3«Y>6«Z) 
1640 8^R-106*8IN(2*(Y-Z))-01*TJ*8I1I(Y)-84*TJ*C08(Y)+60*8IN(2*Y-3*Z) 

1660 RfR-67*8IN(Y-6*Z)+66*8IN(3*Y-6*Z)+63*8IN(Y-Z)-61*C08(2*Y-3*Z) 

1660 RFB-46*8INa)-28*C08(Y-6*Z)+27*C08(Y“2*Z)-22*C08(3*Y)-21*8IN(2*Y-6*Z) 
1670 L>L«80-»X:d-B*S0 

1580 AN(4)-5. 202561 :EC(4)-.0484042510:M(4)-A(10) 

1600 AP(4)»( (7 .0406E-04*TJ+ . 6004316700000001) *TJ+273. 2776417) *RD: 
ZM(4)-1.3031*RD 

1600 AN(4)»( (3 . 52222E-04*T> 1 . 01063) *TJ-'’00 . 44338611 ) *RD : RM(4)»1047 . 355 
1610 GQ8UB 100 :RP-R: RETURN 
1620 PRINT B<(5):REN 8ATURN(5) 

1630 W»A(10) :X=A(21) : Y-A(22) : Z=A(25) 

1640 L»23046*8IN(Y)+6014*TJ-2680*C08(2*W-6*Y)+2607+1177*8IN(2*W-6*Y) 

1660 L»L-826«C08(2«W-4«Y)-^802*8IN(2«Y)-^425*8IN(W-2*Y)-220«TJ*C08(Y) 

1660 L»L-163*C08(2*W-6*Y)-142*TJ*8IN(Y)-114*C08(Y)+101*TJ*8IN(2*W-6*Y) 

1670 L«L-70*C08(2*X)+67*8IN(2*X)+66*8IN(2*W-6*Y)+60*TJ*C08(2*W-6*Y) 

1680 >L+41*8IN(W-3*Y)+30*8IN(3*Y)+31*8IN(W-Y)+31*SIN(2*(W-Y)) 

1600 L«L-20*C08(2*W-3*Y)-28*8IN(2*W-6*Y+3*Z)+28*C08(W-3*Y) 

1700 L«L+22*TJ*8IN(2*W-4*Y)-22*8IN(Y-3*Z)+20*8IN(2*W-3*Y) 

1710 L«L+20*C08 (4 *W- 10*Y) + 10*C08 (2* Y-3*Z) ♦ 10*8IN (4* W- 10*Y) 

1720 L»L-17*TJ*C08(2*Y)-16*C08(Y-3*Z)-12*8IN(2*W-4*Y)+12*C08(W) 

1730 L=L-12*8IN(2*(Y-Z))-11*TJ*8IN(2*Y)-11*C08(2*W-7*Y) 

1740 L*L+10*8IN(2*Y-3*Z)+10*C08(2*(W-Y))+0*8IN(4*W-0*Y) 

1760 L«L-8*8IN(Y-2*Z)-8*C08(2*X+Y)+8*C08(2*X-Y)+8*C08(Y-Z) 

1760 L»L-8*8IN(2*X-Y)+7*8IN(2*X+Y)-7*C08(W-2*Y)-7*C08(2*Y) 

1770 L=L-6*TJ*aiN(4*W-10*Y)+6*TJ*C08(4*W-10*Y)+6*TJ*SIN(2*W-6*Y) 

1780 L-L-6*aiN(3*W-7*Y) -6*C08 (3* (N-Y) ) -6*008(2* (Y-Z) ) +6*8IN (3*W-4* Y) 

1700 L*L+6*8IN(2*W-7*Y)+4*8IN(3* (W-Y) )*4*8IN(3*W-6*Y) 

1800 L-L+4*TJ*C08(W-2*Y)+3*TJ*C08(2*W-4*Y)+3*C08(2*W-6*Y*3*Z) 

1810 L-L-3*TJ*8IN(2*X)+3*TJ*C08(2*W-e*Y)-3*TJ*C08(2*X) 

1820 L*L+3*C08(3*W-7*Y)+3*C08(4*W-0*Y)+3*8IN(3*W-6*Y) 

1830 L=L+3*8IN(2*W-Y)+3*8IN(W-4*Y)+2*C08(3*(Y-Z)) 

1840 L»L+2*TJ*SIN(W-2*Y)+2*8IN(4*Y)-2*C08(3*W-4*Y) 

I860 LaL-2*C08(2*W-Y)-2*8IN(2*W-7*Y+3*Z)+2*C0S(W-4*Y)+2*C08(4*W-ll*Y) 

1860 L-L-2*8IN(Y-Z) 

1870 B=8207*8IN(Y)-3346*COa(Y)+462*SIN(2*Y)-180*C08(2*Y)+186 
1880 B»B+70*TJ*C08(Y)-71*C08(2*W-4*Y)+46*8IN(2*W-6*Y) 

1800 B-B-46*C08(2*W-6*Y)+20*8IN(3*Y)-20*C0a(2*W-3*Y) 

1000 B-B+18*TJ*8IN(Y)-14*C08(2*W-6*Y)-11*C08(3*Y)-10*TJ 
1010 B»B+0*8IN(W-3*Y)+8*8IN(W-Y)-6*8IN(2*W-3*Y)+6*8IN(2*W-7*Y) 

1020 >B-6*C08(2*W-7*Y)+4*aiN(2*W-6*Y) -4*TJ*8IN(2*Y) 

1030 B=B-3*C08(W-Y)+3*C08(W-3*Y)+3*TJ*8IN(2*W-4*Y) 

1040 >B+3*8IN(¥-2*Y)+2*8IN(4*Y)-2*C08(2*(W-Y)) 

1060 8^66774-B3262*C0S(Y)-1878*8IN(2*W-4*Y)-1482*C08(2*Y) 

1060 R-R+817*SIN(W-Y)-630*C08(1-2*Y)-624*TJ*8IN(Y) 
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1980 

1900 

2000 

2010 

2020 

2030 

2040 

2050 

2060 

2070 

2080 

2090 

2110 

2130 

2160 

2170 
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2230 

2240 

2260 
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2270 

2280 

2290 

2300 

2310 

2320 

2330 

2340 

2360 

2360 

2370 

2380 

2390 

2400 

2410 

2420 

2430 

2440 

2450 

2460 

2470 

2480 

2490 

2600 

2510 

2620 

2630 



BFR-i>349*SIN(2*W-6*Y)i>347*SIN(2*W-6*Y)i>328*TJ*C0S(y) 

RFa-226*8HI(Y) +149*C08(2*W-6*Y) - 126*C08(2* (W-Y) ) 

>8+104*008 (W-Y)+101*C08(2*W-6*Y)+98*C08(W-3*Y) 
B^8-73*C08(2*W-3*Y)-62*C08(3*Y)+42*8IN(2*Y-3*Z) 

>R+41*8IN(2* (W-Y) ) -40*8IN(W-3*Y)+40*C08(2*W-4*Y) 

a^8-28*TJ-23*8IN(W)+20*8IN(2*W-7*Y) 
f af • t QbT • 

AM(6)-9. 664747000000M1 :EC(6)-. 06654609270000001 :N(5)-A(22) 

AP(5)a( (9 . 78642E-04*TJ+1 . 08622069)*TJ+338 .3077722) *RD : IN(5)»2 . 4886*RD 
AN(6)-( (- 1 . 52181E-04*TJ+ . 8731961400000001) *TJ+1 12 . 7903889) *R0 : 
R)l(6)-3498.6 
G08QB 190:RP-R:RETDRN 

CLS : INPUT "Ym(4 DIGITS) ";K 

INPUT “MONTH NUMBER “;M 

INPUT "DAY OF THE MONTH ";I 

INPUT "ZULU TIME (HH.MM88) ";Y$:GOSUB 420:UT=X 

INPUT "UT (DD.MMSS) (+N/-S) ":V3:G0SUB 420:L>X 

INPUT "LON (DOD.NMSS) (+W/-E) ";V$:GOSUB 420:LG»X 

INPUT "TEMP (DEG CELSIUS) ";DC 

INPUT "PRESSURE (MILLIBARS) ";MB 

GOSOB 2720 

JD-367*K-INT(7*(X+INT((M+9)/12))/4)+INT(275*M/9)+I+1721013.6 

CLS 

TXJD-2416020)/36626:DT-81.94*TJ-16:TS-JD-2451646+DT/3600/24 
T0=T8/36626 : TS=T8+UT/24 : TJ=T8/36626+l 
TM|-"(UT):":IF DT-OTHEN TM$*"(TDT):" 

G08UB 490 

Ca(-(2 . 58622E-06*T0+2400 . 061336) *T0+6 . 697374660000001+1 . 002737909*UT 
M0»24 : Qi»FNM((ai) : GS=<a(- . 00029*SINa(5) ) 

DL-(- 17 . 23- . 02*TJ) *8IN(A(6) ) - 1 . 27*8IN(2*A(7) ) + . 21 *8IN (2*A(6) ) 
DL=DL-.2*SIN(2*A(1)):0L=DL/3600*RD:BEM DL=NUTATION OF LONGITUDE 
QE«(( . 00181+TJ- . 0059) *TJ-46. 846) +TJ+84428 . 26+9 . 21*C0S(A(5) ) 

OE=(OE+ . 562*008 (2*A(7) ) ) *8D 

FOR N«OTO 6:0N N+IGOSUB 620.730.1110.1220.1360,1020 
IF N=OOR NalTHEN X=XH:Y=YH:Z=ZH:IF N=ITHEN U^:V^:W=0:GOTO 2410 
IF N-OTHEN U-UH:V-VH:W^WH:XS»X:YS-Y:Z8-Z:U8-U:V8-V;W8-W:G0T0 2410 
X-X8+XH : Y^8+YH : Z-ZS+ZH : U-US+UH : V-V8+VH : W-WS+WH 
A»-DL:GOSUB 280 

PRINT DN$:"TRUE POS, TRUE ECL * EQNX OF DATE ";TN|:G08UB 320:PRINT L;B 
GD=R:L1=L:B1»B:IF N^OTHEN LO>L 
>173.142:BC-R/C:RE31 VELOCITY OF LIGHT 
X=X-8C*U ; Y=Y-RC*V : Z=Z-8C*W 

PRINT :PRINT "APP POS. TRUE ECL k EQNX OF DATE ";TM|: 

G08UB 320 .PRINT L;B 

A=-0E:GO8UB 290:G08UB 320:LP-L:BP«B:IF N=OTHEN LS»L:B8»B 
L-L/16 

PRINT :PRINT "APP POS, TRUE EQU k EQNX OF DATE ";TM|:PRINT L;B 
MO-360 : GH-FNM( 15* (GS-L) ) 

PRINT : PRINT "GREENWICH HOUR ANGLE k DBCL ■;TM| 

X»GH:G08UB 360:L|-V$:X-B:G08UB 360:PRINT L$+"'";SPC(3) : V$+"^" 
LH«(GH-LG)*RD:CL^08(LH) :SL»4IN(LH) :BD-B*RD:SB-8IN(BD) :CB-C0S(BD) 
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2640 U>LT*RD:a-C0S(U»:8>SIlI(L0):X-6X*SB-»CX«CB*CL:60SUB 160: 
Ad-A8:AL«A8/BO 

2660 >-CB*8L:X-8B*a-CB*SX«CL:G08UB 120:AZ-A2/RD 

2600 PRINT :PB1NT 'ALTITODE » ■;AL: :XaAL:6aSUB 360:PRINT 8PC(5) 

2606 PRINT 'AZIMUTH - ";AZ; :X»AZ:G080B 3«0:PRINT SPC(6) 

2670 ON NG08UB 2780,2780.2780,2860.2870 

2680 IF N^ITHEN PRINT 'PHASE* ;X;' AGE a';A6;' DAYS 

2600 IF N>1THEN N0»10: PRINT 'MA(aHTUDE » ';FNR(NG) 

2600 IF 1W)THEN 8«16.004/R:M0»10:PRINT *8D -*;8IR$(FNR(8))+*'* 

2610 IF NOITHEN 2660 

2620 D»23484.8*R:8=036.75/D:M0»10:PRINT *8D »';8TR$(FNR(S) )+*'*; 

2630 8«8*(1+8IN(A0)/D):PRINT ' 8D AUG -■;STR$(FNR<3))+*'' 

2640 X-CQ8(A0)/D:GOSUB 160:PA-AS/RD*60 
2660 PRINT *P IN A - ';STR$(FNR(PA))+*'' 

2660 RF-l/TAN((ALt7.31/(ALt44.4))*RD) 

2670 RF-BF* < (MB-80) /030) / ( 1+ . 00008* (RF+30) * (DC- 10) ) 

2680 BF-RF-.06*SIN((14.7*RF*13)*B0):N0-10;RF-FNR(BF) 

2690 PRINT 'REFRACTION - *;3TR$(RF)+*'* 

2700 GUSUB 2720 tCLS :NEXT N:G0TQ 2090 
2710 : 

2720 PRINT : PRINT 'PRESS C TO CONTINUE* 

2730 C$-INlEr$ :IF C$>«"TmW 2730 
2740 IF C$»'C*THEN RETURN 
2760 GOTO 2730 
2770 

2780 OSIN(BS*RD) *SIN(BP*RD) *€0S(B8*RD) *COS(BP*RD) *COS( (LS-LP) *RD) 

2790 X^:GOSDB 160:Y-RR*SIN(AC) :X»GD-BR*CD:GOSUB 120:PA«ABS(A2/RD) 

2800 ON NGOTO 2810,2840,2860 

2810 K-(l*C08(A2))/2:SN-8IN((LO-Ll)*RD):X-8qR(K):GOSUB 150 

2820 AG-29. 53*AS/PI: IF 8N>0THEH AG-29. 53- AG 

2830 MO-lOO: K-FNR(X) :M0-1 :AG-FNR(AG) : RETURN 

2840 MG-(PA*PA*4.247E-07+.01322)*PA+2. 171S*L0G(GD*RP)-4:BETURN 

2860 MG-.01486*PA*-2. 1716*L0G(GD*RP)-1 .3:RETURN 

2860 MG-2.1716*L0G(GD*RP)-8.93:RETURN 

2870 B0-(168 . 1176+1 . 394091*IJ) *RD: II-(28 . 0743- . 0127991*TJ) *RD 

2880 NN-( ( . 242202*TJ+3 . 98599) *TJ+126 . 3629) *RD 

2890 J>( ( . 0171656*TJ- . 4649142) *TJ+6 . 912936) *RD 

2900 S0-(- . 239992+TJ-2 . 731279) *TJ+42 . 92039 

2910 CB-C08(B9) :Y-8IN(II)*8IN(B9)+C0S(II)*CB*SIN(L9-B0) 

2920 X-CB*C0S(L9-B0) :GOSUB 120:UP-A2/RD 

2930 D-BP*RD:CB-C08(D) :8B-8IN(D) :8J-8IN(JJ) :CJ»C08(JJ) 

2940 AG=LP*RD-NN:SN=8INUG) 

2960 Y-8J*8B+CJ*CB*8N:X-CB*C08(AG):G08UB 120:UU-A2/R0 
2970 BB-8J*CB*8N-CJ*8B 

2980 MG-2 . 1715*LOG(GO*RP) -8.68+. 044*AB8(UP+8Q-UU) -2 . 6*AB8 (BB) +1 . 26*BB*BB 
2990 RETURN 
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D* Program Aimotation — ^NAVEPHEM 



Line(8) 



Usage 



30 

40 



50 

60 

70-80 

90 

120 

150 

160 

190-200 

210 

220-260 

280 

290 

320-330 

360-390 

420-460 

490-600 

620-720 

720-1100 

1110-1210 

1220-1350 

1360-1610 

1620-2070 



2090-2240 



P2 = t/2, pi = t, TP = 2t, RD = degree-to-radian conversion. 

EK = Sun’s gravitational constant. MBS = atmospheric pressure m millibar’s 
(1010 is the default value), DCS = temperature in degrees Celsius (10 is the 
default value). 

FNA converts angles in revolutions to degrees between 0 and 2r, FNR rounds 
to the nearest MO. 

FNM is the X modulus MO function. 

BS is an array of planet names. 

GOTO 2090 to commence execution. 

Two argument arctangent function with output interval of 0 to 2x. The -lE-9 
must be replaced by +1E-9 in the PC-1500A. 

The arcsine function can be replaced inlin e by the ASN function, m the PC- 
1500A. 

The arccosine function can be replaced inline by the ACS function in the PO- 
1500A. 

Unscale the radius vector by a factor of l£5. Convert spherical coordinates to 
rectangular coordinates. 

Return if the body is the moon. 

MU is the reduced mass of the body. Compute the approximate heliocentric 
velocity vector of the body. 

Z-axis rotation. 

X-axis rotation. 

Convert rectangular coordinates to spherical coordinates. L = longitude, B = 
latitude and R = distance. 

Convert decimal degrees to degrees, minutes and tenth minute notation. 

Convert DDD.MMSS or HHAIMSS format to decimal. 

Compute the fundamental planetary arguments [Ref. 10]. 

Compute the geocentric spherical position and velocity of the Sun. Convert to 
geocentric rectangular. 

Compute the geocentric spherical and rectangular coordinates of the Moon. 

Compute the heliocentric spherical position and velocity of Venus. Convert to 
heliocentric rectangular. 

Compute the heliocentric spherical position and velocity of Mars. Convert to 
heliocentric rectangular. 

Compute the heliocentric spherical position and velocity of Jupiter. Convert 
to heliocentric rectangular. 

Compute the heliocentric spherical position and velocity of Saturn. Convert to 
heliocentric rectangular. Saturn’s longitude and latitude are saved as L9 and 
B9 on line 2030. 

Input prompts should be designed so that default values or previously entered 
values are displayed. Pressing the return key will input a displayed value. If 
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any valne is changed, then the entire value must be keyed in. This feature 
had to be removed on the TRS-^0 Modd.-4 version since it was difficult to 
implement. 

2250 Computation of the Julian Day Number. 

2270 TJ = number of Julian centuries from noon, January 0, 1900 to midnight of 
input date. DT = AT correction factor. [NOTE: This equation should be 
examined for accuracy yearly — see text). TS = number of Julian days from 
noon, January 1, 2000 to 0000 hours Universal Time of the input date. 

2280 TO = number of Julian centuries from noon, January 1, 2000 to 0000 hours 
Universal Time (UT) of the input date. TS = nunffier of Julian days from 
noon, January 1, 2000 to the input date and time. TJ = number of Julian 
centuries from noon, January 0, 1900 to the input date and time. 

2290 If the AT correction is set to zero, then the positions are referenced to TDT 
(terrestial dynamic time) rather than to UT. 

2300 Compute fundamental arguments. 

2310 GM = Greenwich Mean Sidereal Time. 

2320 GS = Greenwich Apparent Sidereal Time. 

2330-2340 DL = nutation of longitude. 

2350-2360 OE = obliquity of the ecliptic corrected for nutation. 

2370 FOR loop to cycle throng Sun, Moon and planets. N = body number. The 
loop en^ on line 2700. 

2380 X, Y and Z are the geocentric coordinates of the Sun or the Moon. Set the 
velocity components of the Moon equal to zero. 

2390 Save the position (XS, YS & ZS) and velocity (US, VS & WS) components of 
the Sun. 

2400 Compute the geocentric position and velocity components of the Nth planet. 

2410 Corr^ position for nutation. 

2420 Display true position for the true ecliptic and equinox of date. 

2430 Save true distance GD, longitude LI, latitude B1 and solar longitude LO. 

2440-2450 Correct position for aberration (light time). 

2460 Display apparent position for the true ecliptic and equinox of date. 

2470 Convert ecliptic coordinates to equatorial coordinate. Save the right ascension 
LP and declination BP of the body. For the Sun, save right ascension and 
declination as LS and BS. 

2480 Convert right ascension from degrees to hours. 

2490 Display apparent position for the true equator and equinox of date. 

2500 GH = Greenwich Hour angle. 

2510-2520 Display Greenwich Hour angle and decimation. 

2530-2565 Compute and display the altitude and azimuth. The altitude is saved as A9 on 
line 2540. 

2570 Subroutine call for N = 1 through 5. 

2580 Display the Moon’s phase and age. 

2590 Display the planet’s magnitude. 

2600 Compute and display the Sun’s semidiameter [Ref. 2, pg. Bl6j. 
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2610 if the body is not the Moon. 

2620 G)mpnte and display the Moon’s semidlameter [Ref. 2, pg. B16]. 

2630 Compute and display the Moon’s augmented semidiameter [Ref. 2, pg. B16]. 
2640-2650 Compute and display the Moon’s parallax in altitude [Ref. 2, pg. Bl6]. 
2660-2690 Compute and display the refraction correction at the body’s altitude [Ref. 2, 
pp. B14-B15]. 

2700 Continue prompt. End of FOE-NEXT loop on body numb^ N. 

2720-2760 Continue prompt query. 

2780-2790 Compute the phase angle PA. CD is the cosine of the elongation [Ref. 4, 
pg. 312]. 

2800 Transfer for the Moon, Venus or Mars. 

2810 Compute the Moon’s phase angle K [Ref. 4, pg. 311]. 

2820 AG is an approximation to the Moon’s age develop^ by the author. 

2830 Round the Moon’s age to the nearest day and return. 

2840 MG is the magnitude of Venus [Ref. 4, pg. 314]. 

2850 MG is the magnitude of Mars [Ref. 4, pg. 314). 

2860 MG is the magnitude of Jupiter [Ref. 4, pg. 314]. 

2870-2980 Compute the magnitude MG of Saturn. This computation is complicated by 
the varying aspect of Saturn’s rings. Details follow. 

2870-2900 Computation of the ring’s ascending node BO, inclination II, right ascension 
NN, inclination to the mean equator JJ and arc SO from NN to BO. These 
equations result from curve fitting the table in Reference 4, page 365. 
2910-2920 Computation of UP = XT' [Ref. 4, pg. 364]. 

2930-2970 Computation of UU = U and BB = sin B [Ref. 4, pg. 345]. 

2980 MG is the magnitude of Saturn [Ref. 4, pg. 314]. 
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APPENDIX: The QATN Pimction 



This routine is the standard arctangent function corrected for quadrant. The quadrant arc- 
tangent function is occasionally implemented as the ATAN2 function, the ANGLE function 
or the Rectangular-to-Polar function. 

Entering variables are the x- and y-coordinates, X and Y. The exiting variable 
the angle 0, where —r < 0 < ir. Use of the quadrant arctangent function is denoted 
0 = qatn(y,X). 

1. If X ^ 0, go to step 4. 

2. Set 0 = (ir/2) * 8gn(y). 

3. Go to step 8. 

4. Set 0 = arctan(y/X). 

5. If X > 0, go to step 8. 

6. Set 0 = 0 + S’ ♦ 8gn(y). 

7. If y = 0, 8et 0 = ir. 

8. Return. 

Note; 

If y > 0 then 8gn(y) = 4-1. 

If y = 0 then 8gn(y) = 0. 

If y < 0 then 8gn(y) = -1. 

Users of Microsoft BASIC can simplfy the qatn function significantly by using the code 
given below. To return an angle of 0 (designated by A in the code) in the range of (-ir, r), 
use: 

PI - 4*ATN(1); TP - PI 4- PI: EPS - lE-33 
A - ATN(Y/(X-EPS*(X-0))) - PI*(X<0)*(SGN(Y) - (Y=0)) 

To return a value of A in the range of (0, 2r), use: 

PI - 4*ATN(1); TP - PI 4* PI: EPS - lE-33 
A - ATN(Y/(X-EPS*(X=0))) - PI*(X<0) 4* TP*(X >= 0)*(Y<0) 
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